{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RT1-P1 Element for Poisson Equation in 2D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example is to show the rate of convergence of the RT1-P1 mixed finite element approximation of the Poisson equation on the unit square:\n",
    "\n",
    "$$- \\nabla (d \\nabla) u = f \\; \\hbox{in } (0,1)^2$$\n",
    "\n",
    "for the following boundary conditions\n",
    "- pure Dirichlet boundary condition: $u = g_D \\text{ on } \\partial \\Omega$.\n",
    "- Pure Neumann boundary condition: $d\\nabla u\\cdot n=g_N \\text{ on } \\partial \\Omega$.\n",
    "- mixed boundary condition: $u=g_D \\text{ on }\\Gamma_D, \\nabla u\\cdot n=g_N \\text{ on }\\Gamma_N.$\n",
    "\n",
    "Find $(\\sigma , u)$ in $H_{g_N,\\Gamma_N}(div,\\Omega)\\times L^2(\\Omega)$ s.t. \n",
    "\n",
    "$$ (d^{-1}\\sigma,\\tau) + (div \\tau, u)  = \\langle \\tau \\cdot n, g_D \\rangle_{\\Gamma_D} \\quad \\forall \\tau \\in H_{0,\\Gamma_N}(div,\\Omega)$$\n",
    "\n",
    "$$ (div \\sigma, v)                =  -(f,v) \\quad \\forall v \\in L^2(\\Omega) $$\n",
    " \n",
    " where \n",
    " \n",
    " $$H_{g,\\Gamma}(div,\\Omega) = \\{\\sigma \\in H(div,\\Omega); \\sigma \\cdot n = g  \\text{ on } \\Gamma \\subset \\partial\\Omega \\}.$$\n",
    "\n",
    " The unknown $\\sigma = d\\nabla u$ is approximated using the incomplete P2 Raviart-Thomas element (RT1) and $u$ by piecewise constant element (P1). The expected convergence rate is 2nd order in H(div) and L2 norm for $\\sigma$ and L2 norm of $u$. \n",
    "\n",
    "**References**\n",
    "See page 12, Table 9.2. Arnold, Douglas N. and Falk, Richard S. and Winther, Ragnar. Geometric decompositions and local bases for spaces of finite element differential forms. Comput. Methods Appl. Mech. Engrg. 198():1660--1672, 2009.\n",
    "\n",
    "**Subroutines**:\n",
    "\n",
    "    - PoissonRT1\n",
    "    - squarePoissonRT1\n",
    "    - mfemPoisson\n",
    "    - PoissonRT1mfemrate\n",
    "    \n",
    "The method is implemented in `PoissonRT1` subroutine and tested in `squarePoissonRT1`. Together with other elements (RT0, BDM1), `mfemPoisson` provides a concise interface to solve Poisson equation in mixed formulation. The RT0-P0 element is tested in `PoissonRT1mfemrate`. This doc is based on `PoissonRT1mfemrate`.    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RT1 Incomplete Quadratic H(div) Element in 2D\n",
    "\n",
    "We explain degree of freedoms and basis functions for RT1 element on a triangle. \n",
    "\n",
    "### Asecond orientation\n",
    "The dofs and basis depends on the orientation of the mesh. We shall use the asecond orientation, i.e., `elem(t,1)< elem(t,2)< elem(t,3)` not the positive orientation. Given an `elem`, the asecond orientation can be constructed by \n",
    "\n",
    "        [elem,bdFlag] = sortelem(elem,bdFlag);  % ascend ordering\n",
    "        \n",
    "Note that `bdFlag` should be sorted as well. \n",
    "\n",
    "The local edge is also asecond `[2 3; 1 3; 1 2]` so that the local orientation is consistent with the global one and thus no need to deal with the sign difference when the positive oritentation is used. Read [Simplicial complex in two dimensions](../mesh/scdoc.html) for more discussion of indexing, ordering and orientation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA\nB3RJTUUH4wcPBBkNVP9yswAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAxNS1KdWwtMjAxOSAxMjoyNToxM9G1IdUAABu+\nSURBVHic7d19dBTlvcDxJyGBSCCsBKJwyckEIaEF2tW0hcRKNsKhRXlrrRyl9mQDVW69alE8p9V6\nDbmn19oe06ZytZVqs+nx+pLWIzW5viDY8aVJECMrISIoZC7LDSSwcQ3GACHJ/WPsmia7SxJ2dp6d\n+X7+2sxsJr/MCfkyO7OThP7+fgEAgNkSzR4AAAAhCBIAQBIECQAgBYIEAJACQQIASIEgAQCkQJAA\nAFIgSAAAKRAkAIAUCBIAQAoECQAgBYIEAJACQQIASIEgAQCkQJAAAFIgSAAAKRAkAIAUCBIAQAoE\nCQAgBYIEAJACQQIASIEgAQCkQJAAAFIgSAAAKRAkAIAUCBIAQAoECQAgBYIEAJACQQIASIEgAQCk\nQJAAAFIgSAAAKRAkAIAUCBIAQAoECQAgBYIEAJACQQIASIEgAQCkQJAAAFIgSAAAKRAkAIAUCBIA\nQAoECQAgBYIEAJACQQIASIEgAQCkQJAAAFIgSAAAKRAkAIAUCBIAQAoECQAgBYIEAJACQQIASIEg\nAQCkQJAAAFIgSAAAKRAkAIAUCBIAQAoECQAgBYIEAJACQQIASIEgAQCkQJAAAFIgSAAAKRAkAIAU\nCBIAQAoECQAgBYIEAJACQQIASIEgAQCkQJAAAFIgSAAAKRAkAIAUCBIAQApJZg8ASOrcuXMej6eh\noeHgwYOZmZmLFi1at25dcnKy2XMBlpXQ399v9gyAdNrb26+77rq33npr4MIvfelLb7zxxpQpU8ya\nCrA2XrIDQvjJT36i1yg5OTkzM1NfuH///vXr15s6F2BlHCEBgx09ejQrK6uvr2/SpEkNDQ1z5syp\nqalZuXKlECI5ObmzszMlJcXsGQEL4ggJGGz//v19fX1CiBtuuGHOnDlCiBUrVlx22WVCiJ6enkOH\nDpk8H2BRXNQAu9M0raqqSlVVIYSiKKWlpWfPnnU6nUKIRYsW6c/p6uo6efKkECIxMTE7O9u8YQEr\nI0iwNY/HU1JSMnCJqqput3vPnj3BJd3d3TfddNMnn3wihLj66qvHjx8f6ykBe+AcEuxLVdWioqKh\nyxVFqaysdLlcQog9e/YUFxc3NTUJIaZOndrQ0DBz5swYzwnYBOeQYF9VVVUhl2uaVlZW1tfX96tf\n/WrBggV6jZxO565du6gRYByOkGBf2dnZmqaFXJWVlbVgwYLq6mohxJgxY+67776f/exnvCsWMBRB\ngn0lJCSc9zkXX3xxTU3NlVdeGYN5AJsjSLCvoqIi/eK6CAoLCxcsWDBwya233pqVlWXgWIBdESTY\n19BL7IbjlVdeWbp0qRHzADbHRQ2wL0VRFEUZujwjIyPmswDgCAl2pWlaSUmJY3VKwPeJWl6vL1QU\nxe12l5aWmjsbYE8ECXak10i5LUMpyBRCBHyd3up9ij+3srLS7NEA++IlO9jOoBoJIRyZaUKIkC/f\nAYgZggR7GVojAJIgSLARagTIjCDBLqgRIDmCBFugRoD8CBKsjxoBcYEgweKoERAvCBKsjBoBcYQg\nwbKoERBfCBKsiRoBcYcgwYKoERCPCBKshhoBcYogwVKoERC/CBKsgxoBcY0gwSKoERDvCBKsgBoB\nFkCQEPeoEWANBAnxjRoBlkGQEMeoEWAlBAnxihoBFkOQEJeoEWA9BAnxhxoBlkSQEGeoEWBVBAnx\nhBoBFkaQEDeoEWBtBAnxgRoBlkeQEAeoEWAHBAmyo0aATRAkSI0aAfZBkCAvagTYCkGCpKgRYDcE\nCTKiRoANESRIhxoB9kSQIBdqBNgWQYJEqBFgZwQJsqBGgM0RJEiBGgEgSDAfNQIgCBJMR40A6AgS\nzESNAAQRJJiGGgEYiCDBHNQIwCAECSagRgCGIkiINWoEICSChJiiRgDCIUiIHWoEIAKChBihRgAi\nI0iIBWoE4LwIEgxHjQAMB0GCsagRgGEiSDAQNQIwfAQJRqFGAEaEIMEQ1AjASBEkRB81AjAKBAlR\nRo0AjA5BQjRRIwCjRpAQNdQIwIUgSIgOagTgAhEkRAE1AnDhCBIuFDUCEBUECReEGgGIFoKE0aNG\nAKKIIGGUqBGA6CJIGA1qBCDqCBJGjBoBMAJBwshQIwAGIUgYAWoEwDgECcNFjQAYiiBhWKgRAKMR\nJJwfNQIQAwQJ50GNAMQGQUIk1AhAzBAkhEWNAMQSQbIXVVU9Ho+qqud9pt1qFPB1qqo6nD0DwCAJ\n/f39Zs8AY2maVlVVFfxtm+6Y5g8cUxRFURSXy1VYWOhyuYZ+ih1qFPB1eqv3aXVHtXqfGN6eAWAc\ngmRlZWVlmzdvFkKkO6alO6bnKHkFzuXpjun+QOtB7d0D2jv13lohhKIobre7tLRU/yw71Egtr1PL\n68UI9wwAQxEkyyopKanZ9kq+c0Wukpej5IV8jj/QKoSo89bWe2tu37ihtLTUDjXatvFl3/aOke6Z\n2M4I2BFBsia9Ru7Vm8P9wh3EH2gt92y4feMGVVXtUKNR7BmaBBiNIFnQSGuk03/zZi6dvLri28bN\nZq6R1khHk4DYGKOfY4BljK5GQojxKROdc1wvP/fn7tOnLHmENLoaiX/smYrfb+4+8ymXOQDG4bJv\nS/F4PKOrkS7dMX2T+7GmJzRvdXPUZzOXt7p5dDXS6XtmS8VjHo8n2qMB+Bwv2VlKdna2y+nOdy6/\nkI3UqFvfOvTMxl03R2sqGVQs+MPqvLsufM/s07a3tLREayoAA3GEZB0ej0fTtBzlishPe7LmgV8+\nsa75o7pwT8hV8gK+Tq3OF+0BTeOtbg74OkPumZ6eM8/v+K/7t3z3tp8X3L/luy++8UTPubPhtpOr\n5GmaxptnAYMQJKtJd0yPsLb5o7q33t122Lf3VFcg/BamCSH094rGWGpH17wXm1I7uozYeMg9s/XP\nP335LU+b/0jPubNt/iN/fe13j//l3vBbmCaEeP31140YD0CS2QMgaqqqqiK8JPU/rz9+yLf3/UMN\n/f19kbejv1FUqzsqNkV7xFBSO7qyGw4LIea/1BRcuO+a+VH8Et5nm0PumUO+vXsPvimEuCrvOwWX\nr9xR/9+NzTu8H6j+wDG9PYPoe0ZVVS63A4xAkKxDVVX36s1h1+6u7vy0Y5ibKnCu2Nb46+iMFca8\nF5vEP0coaOauw+2zM6L4tTI+arvGpTi6Bn+t44dfFUIkJyXfctWyscki+ap1jc07hBCBUydCBkkI\nUeBcoXo9UZwNQBBBsgj96q8IJ5CuLby5+/Sp7jNdr7xVJYRIPtcZYWs5yhWBbZ1anc+g678XP7wj\n48P2cGtTO7oWP7wzml9uUbroe04cfm7Q8hnj+65fMnXcGHGpdv+x7t7X2hYKIcYmp2RNmxNuUznK\nFZ5tm1VV5fpvIOo4h2QpEU4gub5+/bKr1i1Z+H39w6wT1SlnwyYh3TE93THNoNNI815silCjWJqS\nkpgzKSlrQtLThz/71isn6rw1QoiS7/xHUtLYcJ+i7xlOIwFGIEgWoSiKEOKg1hj5aWmf7dcfJPd2\nOlvuidAkf+CYI3NS9Ab8QsvCmV2TU43Y8qhNTUmcnZaUmJgohHh+x5aOT45HeLI/cCwrKytWowE2\nwkt2FqEH6YDWGPmNn5P+ESQhRMrZdmfLPd7sX5weO/iEjX6va+eaudEfVIiuyak7f7xk4ZP1EY6T\nonsOSas7qt/Ve9DyE5+d/qynN21c8pVzxl45Rzw/du2Dj7vbO3y7m1751jeLQ25K3zNutzuK4wHQ\nESSL0P+Ez3mPkHxTvivE74MfhmvSAe0dg2qk65qcuvOOJakdXan+TzM+bL/ko7aBcWpZMLPhpoVR\n/HKe16vH7h23yf3AoOUVf7p1/+G3vzZv6c3fe0AIkdXXOzY55WzP6SPHD4Tb1AHtHWoEGISX7KzD\n5XLpfzRhRPQmDXrt7qDWGIPb2XVNTm2ffcm+a+bvvGPJC2Wrdt6xuGnZ/PbZGW1RPTwSQigFM0Lu\nmUunZgsh3vtAPXLsg76+PnX3n8/2nBZC/EvGrHCbOqg1FhYWRnc8ADqCZB3FxcX+wLHzHiTpTqTl\nBx8PalK9t9YfOKbkx/T+qgPj1LJgZnQ37lwzL+SeWbzgxuTkcT3nzj6w9Qd3/tL17EsPCSEmjHfk\nf/XakNvR9wzX1wEGIUjWETyNFOE5CQmfPziZlq9l3BhcPvQ4yZGZZsSQptC/l6F7ZurkGXd8f8v0\njMv6+/tPn/lMCDE76/KNP3jk4kmXRNiavp8BRB03V7UUj8dz9533bnI/FvkGQkFK21NK+9PBD0+P\nzdh58e13PfrD1RXfNvQcUux5q5vf3PxeuD3T1d35cWdb+qRpF6VMCLcFf6D13oqVlZWVnEMCDMLf\nQ7IUp9PZfebTit9vds5xjU+ZeN7nBybMF0I4uvbpHyb1dqWe2HHyqsvy7ojmNQUyuHRuRvfpU9sq\nnwy5Z8Ymj0ubkJ4c/u1Heo3cbjf/XgDjECSrcblc+w/s8zxTMbompSUnfjM14ehXM3suCvvbOU4p\nBZkRmhRBsEaVlZXGjQeAc0hWo2mapmmZSyeXezYM86I77ZK1jWNdwQ9TO7oW/3aHQbfcNpdrU8H8\n9crw94ygRkAMcYRkKZqmlZSUKLdlLLw5Tz8a+Oz0pwkRbynkD7TubHj6F7VVSkFm7unT+sKx3T0z\n9h619nHSMPfM7565mxoBscFFDdYRrFHwLURqeZ1aXi+ESHdMy1HycpWv6Q/8gdY6b60/0Krfd8CR\nmeZcM9e1qWDei00Db7+t31JBttv8RMUw94yiKG63mz82AcQGQbKIoTUKCvg6tXqf99nmQTdL1Tuk\n5GcO/BT7NEmE3zN6hwoLC3nLERBLBMkKItRooICvM+D7RH8c4Zm2apIu4OtUy+uc479eXFxMhACz\ncC+7uDfMGgkhHJlpw3m7q/7XWoNN0q9xsHaT9D2jTFSoEWAirrKLb8Ov0Yjsu2Z+07Iv/oi4ha+7\nAyAPghTHDKqRjiYBiDGCFK8MrZGOJgGIJYIUl2JQIx1NAhAzBCn+xKxGOpoEIDYIUpyJcY10NAlA\nDBCkeGJKjXQ0CYDRCFLcMLFGOpoEwFAEKT6YXiMdTQJgHIIUBySpkY4mATAIQZKdVDXS0SQARiBI\nUpOwRjqaBCDqCJK8pK2RjiYBiC6CJCnJa6SjSQCiiCDJKC5qpKNJAKKFIEknjmqko0kAooIgySXu\naqSjSQAuHEGSSJzWSEeTAFwggiSLuK6RjiYBuBAESQoWqJGOJgEYNYJkPsvUSEeTAIwOQTKZxWqk\no0kARoEgmcmSNdLRJAAjRZBMY+Ea6WgSgBEhSOawfI104ZqU2tG1+OEdGR+2mTgbANkkmT2AHdmk\nRrp918wXQsx/qUn/cOBx0nzRtHP2JWYOB0AmHCHFmq1qpBt6nKQ/yPiwPXvXYZOGAiAdghRTNqyR\nblCTgua/2BT7YQDIiSDFjm1rpGufnTF0YWpH18InG2I/DAAJEaQYsXmNhBDhwpPxYRtXNwAQBCk2\nqFFqR9fhBTPDrQpe8gDAzgiS4aiREKJrcuq+a+a/ULYq5Jkkrm4AIAiS0ajRQHqWdt6xuGty6qBV\nMwkSYHsEyUDUKKT22Zfs/PGSQYdKXZMnmDUPAEnwxlijUKMI9EOlloUzsxsO68dGbaGuwQNgKwTJ\nENRoOPQs6bdyAABesos+agQAo0CQoowaAcDoEKRookYAMGoEKWqoEQBcCIIUHdQIAC4QQYoCagQA\nF44gXShqBABRQZAuCDUCgGghSKNHjQAgigjSKFEjAIgugjQa1AgAoo4gjRg1AgAjEKSRoUYAYBCC\nNALUCACMQ5CGixoBgKEI0rBQIwAwGkE6P2oEADFAkM6DGgFAbBCkSKgRAMQMQQqLGgFALBGk0KgR\nAMQYQQqBGgFA7BGkwagRAJiCIP0TagQAZiFIX6BGAGAigvQ5agQA5iJIQlAjAJAAQaJGACCFJLMH\nMJndaqTV+QJHOx0z0mzy/QKII7YOkk1qFPB1eqv3aXVHtXqfECLdMc0fOObITHPMmKQUzFDyM639\n7QOIF/YNkh1qpJbXqeX1Qoh0x7R0x/TlrmUFzuXpjun+QOtB7d0D2jtqea0Q9Y7MNOeaua5NBWbP\nC8DWbBokO9Ro28aXfds7lrtuyVXycpS8gavSHdPzndPznctXuG4RQtR5a+ufqBFC0CQAJrJRkE6d\nOvXII480NjYeOXLkyJEj89bPtnyN3Ks3D0rRIOmO6UKIFa5bCpzLy5/YIGhSKA8++OCpU6eKioqW\nLFli9iyAldklSO+///6yZcuOHDkSXHL8P4//357j1z+2InFMgomDGWGYNRoo3TF9k/sxmjTU22+/\nfc899wghuru7CRJgKFtc9t3X17du3Tq9RomJiZMVR0KCEELsf/HDhj80mjxctI2iRjq9SU1PaGp5\nnUGzxZdz58698MILq1atMnsQwC5sEaSmpqZdu3YJIVJTU9c8seKOuvWrf7tMX/XRay2mjhZl3urm\n0dVIF2ySt7o56rPFl+Li4gkTJqxater48eNmzwLYhS1esnv33Xf1B7OvVeZ8a5YQYvbV2Z+vS7DU\n63Vqed1q110ha9Qv+l9reGbX3peOnTg8MXXy3Fn5K4s2TEydPOhp6Y7p+c4VavkzzjVzYzKypPx+\n/5kzZ8yeArAXCwZJ07SqqipVVYUQiqKUlpYuWrRo0aJFWT+aOvOqrC5/d0fLx39/dLf+5FkuxcRR\no8tb3RzwdeZcd0XItc++9NDfdj2rP/YHWt9457mPjnjvvflPycnjBj0zV8mrVbdqdT4LX/QxkP4+\nLW91szbmmKZppaWliqLceeed119/vRCirq5u69atZs8I2ILVguTxeEpKSgYu0cvkeugb+q/XZ9zb\nfI2t+qp5K3O/sf7ymM9oLP3CuUFOftz6+u6/CCHmzspfvHBtvbdm977tre2HGt/fsfCr1w7ZwjQh\nhFZviyB5q5u3bXxZfxwQnR6PR1VVt9tdWlqqLxwzZgxBAmLDUueQVFUdVCMhhKZpgd6O4IcJiV+8\nRveRqjW/cCBGwxnP+2xzvnN5yFWHfN6+vt6EhMTiVaVzZ+WvXX7PmDFJQoijbR8OfXK6Y3qOkqfV\nHTV2XAlodb5gjb5YqGl6lsyYCLA1SwWpqqoq5PKAr1O/YYEQYt1fb/jpB7dd859XJySI051nXrrv\ntd6evhjOaCCt3perfC3kqnPneqZOnpE9Y96kiVOEEP19ff39fUKIoeeQdAXOFYGjnxg3qiTCXbuh\naVpZWVmMhwFgqZfsIvyv1tfYuuWbf1xaWjgudawQIiN3yqQZaQFf5+nOM29X7pk2LyN2UxpDv09d\njhL6BNKVV6y68orPL18+23P6j8/f39fXl5CQOPeyhSGfn6NcEdjW6a1udsxIM2hgGeg7LfQqTYvh\nIACEsFiQIvwS6T3b6z/88Xu/PJiRkSGEOHfu3KnWLn3VyZpTZ9SYzGckTTsmwpxACrr0452fHt72\n72/uPRzoEkKsLPrXGZfmhHxmumN6umOat+IDRVEMGFYWAV9nuFUECYg9SwXJ5XJFfum/ra1t06ZN\nF110UWVlZW9vrxAiLS2toaEhKSnu94OqqkVFRQe1xgjvQFLf2/nomw1nevvHJyX8aNHSnEXrImzQ\nHzj20G8q3W539GeVRlFRUbgfGGuXGJCTpc4hFRcXh1uVm5srhPD7/evXr1+7du2rr76qL3/88cct\nUCPxj1+gB7SwN57Y2fDUb9TtZ3r7v+xI3rZ4ytLcSG8zqvfWCiGsXSMR8QcmeJUdgJixVJBcLpfL\n5Qq5fPfu3XfffXdKSkpw4Ve+8pXa2lr9vSYWoCiKy+U6GCZIbSf/9y/bfyuEmDwusezySef6xZGP\nO46daAl0tod8/gHtHcvXSET8gQl++4mJlvo3Asgsob+/3+wZokl/V+zmzZv1DxVFGfiekt7e3paW\nFr/fP2vWrPT0dNOmNEZZWdmWisce2FgzdNVzrz68/e9/Grr86/OW/vB7Dwxdfm/Fiod+84AdmhT5\nBwZALFktSLrgGWlbnQnQNC07O3uT+7Ghp5Eeffqu9w68MfRTvjH/2+uv+/mghfXeWs+2zS0tLfbZ\ne/b8gQFkY4XTJ0PZ89dK8DTS0CDdeuOvhRBK21NK+9P6Ei3jRu2Stefdmk3Y6psFpMXr45ZSWVlZ\n763xB1pHvQV/oNWzbXNlZWUUpwKA4SBIluJ2u2/fuKHcs2F0TfIHWu+tWOl2u+1w9giAbAiS1ZSW\nlo6uScEacXgEwBQEyYJG0SRqBMB01ryoAfqFy+UVG/KdK3KVvAi3b/AHWuu8tbXqVmoEwFzWvOwb\nurKyMv0dNumOaTlK3tppp5aOP6iv2v5ZzlPHJup3ZODNNwBkQJCsT9M0VVWrqqrmtb1z+5cn6gu3\nvH+qtnuK2+0uLCwMebcCAIgxgmQj/upyf3W5/jh9zab0NZvMnQcABuKiBgCAFAgSAEAKBAkAIAWC\nBACQAkECAEiBIAEApECQAABSIEgAACkQJACAFAgSAEAKBAkAIAWCBACQAkECAEiBIAEApECQAABS\nIEgAACkQJACAFAgSAEAKBAkAIAWCBACQAkECAEiBIAEApECQAABSIEgAACkkmT0AYqGn3depVnf+\nrTq45LPmuuS/ZSZlzBg/t8DEwQAgKKG/v9/sGWA4X+l13c31Q5enudZceltF7OcBgKF4yc4W0tds\nCrmcwyMA8iBItpA8NTPk8ovm5sd4EgAIhyDZQnJGZsj2JGeEDhUAxB5Bsq801xqzRwCALxAkuxh6\nGokTSACkQpDsYuhpJE4gAZAKQbKLoaeROIEEQCoEyaY4gQRANgTJRi79twr9hbvkqZlpRQQJgFy4\nU4O99LT7/NXl3J0BgIQIEgBACrxkBwCQAkECAEiBIAEApECQAABSIEgAACkQJACAFAgSAEAKBAkA\nIAWCBACQAkECAEiBIAEApECQAABSIEgAACkQJACAFAgSAEAKBAkAIAWCBACQAkECAEiBIAEApECQ\nAABSIEgAACkQJACAFAgSAEAKBAkAIAWCBACQAkECAEiBIAEApECQAABSIEgAACkQJACAFAgSAEAK\nBAkAIAWCBACQAkECAEiBIAEApECQAABSIEgAACkQJACAFAgSAEAKBAkAIAWCBACQAkECAEiBIAEA\npECQAABSIEgAACkQJACAFAgSAEAKBAkAIAWCBACQAkECAEiBIAEApECQAABSIEgAACkQJACAFAgS\nAEAKBAkAIAWCBACQAkECAEiBIAEApECQAABSIEgAACkQJACAFAgSAEAKBAkAIAWCBACQAkECAEiB\nIAEApECQAABSIEgAACkQJACAFAgSAEAKBAkAIAWCBACQAkECAEiBIAEApECQAABSIEgAACkQJACA\nFAgSAEAKBAkAIAWCBACQAkECAEiBIAEApECQAABSIEgAACkQJACAFAgSAEAKBAkAIAWCBACQAkEC\nAEiBIAEApECQAABSIEgAACkQJACAFAgSAEAKBAkAIIX/B0sy0uxRZTJ0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "node = [1,0; 1,1; 0,0];\n",
    "elem = [1 2 3];\n",
    "edge = [2 3; 1 3; 1 2];\n",
    "figure;\n",
    "subplot(1,2,1)\n",
    "showmesh(node,elem);\n",
    "findnode(node);\n",
    "findedge(node,edge,'all','rotvec');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Local bases of RT1 element\n",
    "\n",
    "Suppose `[i,j]` is the k-th edge. The two dimensional curl is a rotated graident defined as $\\nabla^{\\bot} f = (-\\partial_y f, \\partial _x f).$ For RT0, the basis of this edge along with its divergence are given by\n",
    "\n",
    "$$ \\phi_k = \\lambda_i \\nabla^{\\bot} \\lambda_j - \\lambda_j \\nabla^{\\bot} \\lambda_i. $$\n",
    "\n",
    "For BDM1, one more basis is added on this edge\n",
    "\n",
    "$$ \\psi_k = \\lambda_i \\nabla^{\\bot} \\lambda_j + \\lambda_j \\nabla^{\\bot} \\lambda_i. $$\n",
    "\n",
    "For RT1, two more bases are added to the triangle. Suppose $i,j,k$ are the vertices of the triangle and $i<j<k$. The two bases are\n",
    "\n",
    "$$ \\chi_1 = \\lambda_j\\phi _{ik} = \\lambda_j(\\lambda_i\\nabla^{\\bot} \\lambda_k -\n",
    "\\lambda_k\\nabla^{\\bot}\\lambda_i),\\quad\n",
    "   \\chi_2 = \\lambda_k\\phi _{ij} = \\lambda_k(\\lambda_i\\nabla^{\\bot} \\lambda_j -\n",
    "\\lambda_j\\nabla^{\\bot}\\lambda_i).$$\n",
    "\n",
    "Inside one triangle, we order the local bases in the following way: \n",
    "\n",
    "$$\\{\\phi_1,~\\,\\phi_2,~\\,\\phi_3,~\\,\\psi_1,~\\,\\psi_2,~\\, \\psi_3, ~\\, \\chi_1, ~\\, \\chi_2\\}.$$\n",
    "\n",
    "Note that $RT_0 \\subset BDM_1 \\subset RT_1$.\n",
    "\n",
    "The first 3 dual bases are the line integral over orientated edges\n",
    "\n",
    "$$d_i^{\\phi}(u) = \\int_{e_i} u \\cdot n_i \\, ds,$$\n",
    "\n",
    "and the second 3 dual bases are\n",
    "\n",
    "$$d_{ij}^{\\psi}(u) = 3 \\int_{e_{i, j}} \\boldsymbol{u} \\cdot \\boldsymbol{n}_{i, j}\\left(\\lambda_{i}-\\lambda_{j}\\right) \\, ds,$$\n",
    "\n",
    "where $n_i = t_i^{\\bot}$ is the rotation of the unit tangential vector of $e_i$ by $90^{\\deg}$ counterclockwise. Note that both the bases and the dual bases depend on the orientation of the edge not necessiliarly couterclockwise.\n",
    "\n",
    "It is straightforward to verify the \"orthogonality\":\n",
    "$$d^{\\phi}(\\psi) = 0, \\quad d^{\\psi}(\\phi) = 0, \\quad d_i^{\\phi}(\\phi_j) = \\delta_{ij}, \\quad  d_i^{\\psi}(\\psi_j) = \\delta_{ij}.$$\n",
    "\n",
    "But it is diffcult to find the \"orthogonal\" dual basis for $\\chi_1,\\chi_2$. One choice is $$l_f^1(u) := \\int_{T} u\\cdot n_{ik}\\, dx, \\quad l_f^2(u) := \\int_{T} u\\cdot n_{ij}\\, dx.$$ Then $l_f(\\phi) \\neq 0, l_f(\\psi) \\neq 0$. The good news is $d_e^{\\phi}(\\chi) = d_e^{\\psi}(\\chi) = 0$ as $\\lambda_k\\phi _{ij}\\cdot n_e = 0$ either $\\lambda_k|e = 0$ or $\\phi_{ij}\\cdot n_e = 0$. \n",
    "\n",
    "We need to solve a small system to get the coefficients for bases $\\chi$. Fortunately this has been done in quadratic H(curl) element. More precisely, suppose $$u=\\sum \\alpha_{i} \\phi_{i}+\\sum \\beta_{j} \\psi_{j}+\\sum \\gamma_{k} x_{k}$$\n",
    "\n",
    "The interpolation is implemented in `faceinterpolate`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Local to global index map\n",
    "\n",
    "Three local edges are `locEdge = [2 3; 1 3; 1 2]`. The pointer from the local to global index can be constructured by\n",
    "\n",
    "    [elem2edge,edge] = dofedge(elem);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "elem2edge =\n",
      "\n",
      "  8�3 uint32 matrix\n",
      "\n",
      "    8    3    2\n",
      "   11    6    5\n",
      "   15   10    9\n",
      "   16   13   12\n",
      "    5    3    1\n",
      "    7    6    4\n",
      "   12   10    8\n",
      "   14   13   11\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA\nB3RJTUUH4wcPBCA64sZaBgAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAxNS1KdWwtMjAxOSAxMjozMjo1OIK8LFAAACAA\nSURBVHic7N15XBNn/jjwT0IgGCBErgokMniggkc8qkKrhNZaa7Xa1mLtJWprj7WtW7f7bbe7K373\nu92fa9nFHmtrPUKPbWV7YLVVW6vBIohnqka8GQwEBYIhQCAQwu+Ph40hF4dJZkI+75evV8PMZHie\nfsjnM5nnmRlOZ2cnIIQQQkzjMt0AhBBCCAALEkIIIZbAgoQQQogVsCAhhBBiBSxICCGEWAELEkII\nIVbAgoQQQogVsCAhhBBiBSxICCGEWAELEkIIIVbAgoQQQogVsCAhhBBiBSxICCGEWAELEkIIIVbA\ngoQQQogVsCAhhBBiBSxICCGEWAELEkIIIVbAgoQQQogVsCAhhBBiBSxICCGEWAELEkIIIVbAgoQQ\nQogVsCAhhBBiBSxICCGEWAELEkIIIVbAgoQQQogVsCAhhBBiBSxIyI+01x3R7Z/RXneE6YYgpzBG\n/ozT2dnJdBsQ8ob2uiONxUsAgCsQDxr1Kl+yiOkWIVsYIz+H35CQv2hTf01emA2VxmtfM9sY5BDG\nyM9hQUL+ol176yyQuaWSwZYgZzBGfg4LEvIL7XVHzIZbCc5sqMRRCrbBGCEsSMgv2B9um7SY7NgF\nY4SwICG/YKortVnSbrcEMQtjhLAgIb/QbnesjUMUbIMxQliQ0MBnMzhB4BAFq2CMEGBBQv7A2YE2\nDlGwB8YIARYk5A/sBycIHKJgD4wRAixIyB/YD04QOETBHhgjBFiQ0IDncHCCwCEKlsAYIQILEhrg\nXB9i4xAFG2CMEIEFCQ1wzgYnCByiYAOMESKwIKEBztngBIFDFGyAMUIEFiQ0kLkYnCBwiIJxGCNk\ngQUJDWS9ObjGIQpmYYyQBRYkNJBZD05wBWKH2+AQBbMwRsgCCxIayKwHJwIjp1teWz+KFIcomIUx\nQhZYkNCAZTM4wRXEO3yNQxQMwhgha1iQ0IBlfVhtfbhN8KwOxnGIgikYI2QNCxIasLiDxHzJIpLU\neFHTbNYGRk0DAK5APGjUq86GLpCnYYyQNR7TDUDIUwKjpgdGTQcAs6GSKxC3XMi1XsuXLBo0ajVD\nTUNdMEbIGn5DQgOfw4NrPOJmFYwRAixICCGEWAILEkIIIVbAgoQQQogVsCAhhBBiBSxICCGEWAEL\nEkIIIVbAgoQQQogVsCAhhBBiBSxICCGEWAELEkIIIVbAgoQQQogVsCAhhBBiBSxICCGEWAELEkII\nIVbAgoQQQogVsCAhhBBiBSxICCGEWAELEkIIIVbAgoQQQogVsCAhhBBiBSxICCGEWAELEkIIIVbA\ngoQQQogVsCAhhBBiBSxICCGEWAELEkIIIVbAgoQQQogVOJ2dnUy3ASEPomk6Ly9PoVBMHaL8nyUi\nsnBHiSj/iEgmk6Wnp8tkMkYbiDBGqAt+Q0ID1rp16zgcTmJi4nu5H2noxlHUZMsqPkQadaHZ2dkZ\nGRmJiYnr1q1jsJ3+DGOErOE3JNRnjY2NH3zwwYkTJ2pqaoYPH/7YY4898MADTDfK1rJly3YV7EuV\nzh9FTU6iJgNA8pD85Nh8svZcdea565lanQYAipW7S5S7Xl79/Nq1a5lssbsdOHDg559/Li4ujoqK\nmjlz5sqVK/l8PtON6sbPY1RdXb158+YTJ050dnbefffds2fPnjhxItONYhgWJNQ3586de+CBB65d\nu2a98JFHHsnPzw8ICGCqVTZIpstamJ1kdcRtn+wsq7Q6TY78+YGU79asWfOPf/zDeklKSsrhw4fD\nw8OZapINP4/R/v37H3/8ca1Wa1nC5XI/+eSTJ598ksFWMQ5P2aE+MJvNy5cvJ9UoODhYKpVyOBwA\n+Oabb3Jzc5luXReHmc61SFHcmqyP3sv9aGCcF1q/fj2pRnFxcStWrJBIJACgUql+//vfM920Ln4e\no4aGhqeffppUo+HDh6empgKA2WxeunRpUVER061jEhYk1AdnzpwpLS0FgLFjx6rV6lOnTuXl5ZFV\ne/bsYbRpXfqR6YgBk+/a2tpINRo8ePCRI0e2bNlSWlo6aNAgAPj222/ZcEYEY5SXl3f9+nUAmDVr\nVllZWXFx8QcffAAAHR0dH330EdOtYxIWJNQHJ0+eJC8WLVoUFRUFAJbRI/JViVlyubx/mY6w5Du5\nXO7upnnPzp07a2pqAODZZ5+VSCQNDQ2BgYHl5eUVFRWnTp1iPEwYIwA4ffo0efHKK68EBgYCwAsv\nvCASiQDgm2++MRqNTDaOUTymG4B8yTPPPLNkyRIACAoKqq2tvXTp0oYNG8iq+++/n9GmAQCsW7du\nvux5h5nuyK8/vPvJZ4KgrlP2hrbPWtr3AsC89GdTRqRZNosUxaVK569bty4rK8srTXa/X3/9lbyI\niIhISUkpKyvr7OwcMWLE+vXrH3nkEWbbBi5jBAAbPt9XXVv9+mLRvZMG2ay6Vn3+ix/+DgBPPPiG\nr8dIp9ORF5bjAy6Xy+PxAMBgMGg0msTERMYaxyj8hoQcoGl63bp1iYmJZEqu5QxJQEBAcHBwcHAw\nl8tdsGDBXXfdVVBQAACLFy9++eWXGW0yyOVymqaTqEkO1+oaa1RXNcfOG8k/1VXNVfXpq+rTjc06\nmy1HUZNpmlYoFB5v8e1xFiPLfJM333zz3Llz5Bzd5cuXH3300a1btzLWXADoKUaqy8W7ik4fO2/U\nNnTYrDK2tWz5+i0SslZjk6/HaNSoUeTFhx9+aDabAeDLL7+sq6sjC2/cuMFIa9kACxKyRdN0YmJi\ndnY2TdPkx+zs7MTERPKjBZd7649n3759+fn53m2mY5GiOIfLqbjkx+6Z/Px8IfknTZIAAIfDiYmU\n2O0hFgAKCws93dTb4SJGDQ0NZBsej7dx48ZTp069+OKLZMkf/vCHtrY2ptpsYR+j7wu3vPvZK+99\nvtrZEFf+3pwbdRVWe/DtGC1fvjwoKAgAvv/++/Hjx99zzz1PPfWU5Y3+fMoOCxKytWzZMvuFNE1n\nZGRYLykqKrp58+Z7773H4XB0Ot0rr7zS3t7urTY6kJeXlyqd52zt6GFTX8m89+1nI95+NuL/VkQY\nWtsAYJ5s5TDxOJstI0VxSdRklh99u4hRdHQ0+fGpp5565ZVXpFLp+++/P3z4cACoqam5dOmSVxva\nnbMYKY7lqy4Xd3aaHb7r5Lmfi04WBPNDLEt8PUbDhw//+OOPBQIBAKhUqoMHD4aGhlIURbYRi8Xe\nbCerYEFC3cjlcmefc5qmFy1aNGXKlClTpty8eRMARCLRqlWrkpOTAUCn0/3yyy/ebKoNhUIxiprS\nmy0/39948dqN4UMnPJj+rMMN0qTzbb4OsorrGFmOr1NSUsgLLpdLYgQA9fX1nm+gU85i9GD6cwvv\nfen+u5far7rZcOPTXX/lcrlLF3S7/MinY6RQKJ555pkzZ87k5uauXLnynXfeOX78OJkJGRAQ4M8F\nCSc1oD4oKSnRaDQA8Nxzz40bNw4AWltbL168SNbu2LGDqZpEcpOzwQlreoP5/z7VAcDcGcs54HjK\nWRI1SV6QvWzZMstBK6u4TsRnzpwhL7Zu3drc3AwAHR0dhw4dIgt37dp14MABDzfQMRcxkt35GADo\nm+r3FeVZL+/sNG/79k+GFv182crhkvHWq3w6Rlu2bCHjeUuWLHn11VcBQKlUlpWVAcDcuXPZdkMN\nb8KChLqpqKhwsTZQwgUNAMB33+8826Tk8rjXz9aQM3UBQQEXQk5fbDzjnXba0Bn04HwAiTh3PfPc\n9cyvfsyta/gsNnpYysg0Z1tGiuIiRbGKCz9TkWw8VqUvVLpYyxllCi7nt+qN58+f3/r95ogEkeb0\nDTKwFBoTcsxcBI3eamh3vYkRcbTiFe6puQCw55dtF+mTwyXj5858trG523c7n47RyJEjs7OzAeDQ\noUOff/45n89fvnw5WWUZ8PNPWJBQN+np6S7WZryepsgppksq21tNF/ZdsV71yHtzU+Ynebh1TtHF\namW+6iJ9oserW46e2QsAk1NmOft6RGh11TMWz5FmprizlW5Cp6rli9TO1kozU8Y9PPqLrJ0mo0l9\nTKM+piHLuTzu49sWiCfFequZtnofI4ujZ/YAgLbh+vqtWR0dJrLwk51/SR6RumTu7303RjNnzpw7\nd+4PP/xw7dq1GTNmWJa/+OKLc+bM8UoDWQrHkFA3MpnM2a3+qVQJlSZZkvdw2otTePxbhzJ3jIl+\n4pOHGaxGACCShAPABfqE683qblY1NNYBwOjEO11sVqLcDQDszHQAQKVJqFTbyYFdq1IlVJpkeDq1\nJG9hvHQIl9f1Ab9jTNSKnUsYrEbQ6xjZ0+lr6Kpz6utdZ4Zr6tU36mjfjZFMJsvIyPj3v/+9YsUK\ny+0fExIS/vjHP37wwQeMX7nMLPyGhGxt3749IyPD5iS4SCKUrUkFAH5o0Ow/pc/6w0ydusFQ3xKR\nKBIMtr2G0ftEEiGVKrnYU7K7VHEKAHi8IPvJddYu0MdZm+kIkUQIJQ4WkhgBwPCZCcNnJrQ1t9de\n1A6mwn0oRtYeve9VQ6uevG4yNOTvzQGAuTOXjxk2vVj5nS/GiKIocnPY8PDwLVu2/POf/7x06dLg\nwYP99kpYG/gNCdmiKOrgwYPZ2dlkuJikudWlz1Fpt474uAGcCEoknhTLhkxHUGli8qgCFy5fUwLA\nYGFMQICrQ7GL9AnrzrJNweq9ALC69DnZmlSRRAhOYgQAQSGB8ROH+FaMrI1Lunva+Lnk35SU+8jC\n5OHTk6hJPhcjiqKys7PLy8utz0CEhYVNmjQJq5EFFiTkADmOKy8vh64P1a3x/6lV56ZWnXvk/CGh\n0cBcAx2QZo7V6qpdH4BfrTwNAKKwGBfblCh3a3XVzs63MI5kuoW5c0QSoWxN2urS52CgxMj1ySoO\n51ay8sUYlZeXD4wHZ3gUnrJDfTNNU0ZePHyhsCwy4Wh8MrPtsSDHoRdcjpmvfam3t5Mge2MbS6Zz\nvZmPxigsJOKj7OPO3isMvbW27qYGfDxGyCH8hoT6wPqIW2g0TNOULT29Z2rVOQabZG1h7pwS5a4+\nnRSyodVp5AXZ7Mwmvcx0GCMGYTW6TViQUN/o+QLrHy0pjw1nh6SZKeNWUDny5/uX77Q6zR9yH5Jm\nprBwtLxPmQ5jxAisRrcPCxLqAz1fcDTOwfkfodHw8IVCNhyGy9ak9S/fWTIdCxNKX6sRxsj7sBq5\nBRYk1DdlUQlVYdH2y9lzdqgf+W6AZTqMkZdhNXIXLEioz/YnOr2HKUvODlny3S7FZtfz7rQ6zS7F\n5oGX6TBGXoPVyI1wlh3qMz1fUBo3xjKVyx45O8Ts/C4yDXp3zubdAJGi2CRq8ihqCnmh1WmKlbu1\nOg252p9cxGM9bZolbifTYYy8A6uRe3GcPRELIQDgcDjZmjX2y0k66/EQW88XMD7tWKfW0yVq5Q4V\nXdLt3mIiiVCamULutcNU21zofabLjsvBGDGiTzHCTNsbWJCQK84KEgCMqauYVe70qhFrer7g21Hp\nNlO/vE+n1uvUXU9TZWeCs+jTcbezggQYI0/qa4ww0/YGnrJD/VQWlTCmriK+sbbHLdlwdggARBIh\nOy+ltOHGs0AYIw/BM3UegpMaUP+5GDm3wZ75XSzn9kyHMXI7rEaegwUJ9R8ZOXexgc3kY5LyxtS5\negagP/NEpsMYuRdWI4/CgoRuy/koysXAQ1hb8/7EKdYbVIVFl0UleKVpPsZzmQ5j5C5YjTwNCxK6\nLc7uC0AIjYapmnPfjkq3HKSXxrs6WvdbHs10GCO3wGrkBViQ0O2yvy+A9Y9ktPxofHLe+AdK48Y4\nvIOAn/NCpsMY3SasRt6BBQm5gf3IuU2+W3p6j54vYM9zENjDa5kOY9RvWI28BgsScgObkfOwtmab\nQ3KS75hoGqt5M9NhjPoHq5E3YUFC7mE9ci40GsbUVexPnIL5zgXvZzqMUV9hNfIyLEjIPWxGzuMb\na+MbazHfOcNIpsMY9QlWI+/DgoTcxuYU0FTNOT1fgPnOHoOZDmPUS1iNGIEFCbmTZeSc3LKTvMB8\nZ43xTIcx6hHjMfJbWJCQO5GR89K4MXnjH7DM18J8Z8GGTIcxco0NMfJbWJCQmx2NT7afOoz5DtiU\n6TBGzrAnRv4JCxLyEj/Pdz6R6TBGwPoYDWxYkJD3+G2+86FMhzFiuiF+DQsS8io/zHc+l+kwRogp\nWJCQt/lVvvPRTIcxQozAgoQY0L98p1PrPdwuN/PpTIcxQt6HjzBHzCD5blb5ccsDtkm+yxv/gPVm\nOrVemX9Wp9Yr81UAIJIIqVQJlSYRiYVUmoSBdvfaAMh0GCPkZZzOzk6m24DYi8PhZGvWeG7/QqPB\nOt8BgJ4vIPlOkVOsyCkhCykKsrIgPR1oGvLyQKEAABBJhNLMFNmaNM81r9+8memy43IwRv3g5Rhh\npu0NLEjIFU8XJHCS7zLWNUK1miQ4mcz2LTQNNA2FhSCXg2hqCtuOcL183O3pggQYo9uGBamXsCAh\nV7xQkMBRvrtWY5I+W9njG2kaMjLYle+8fxbICwUJMEa3BwtSL+GkBsQ8+/HzoTE83X5xj2+kKDh4\nEHRHVSTFMG4Aj0lgjJAXYEFCrEDynUoXZFliNvhYvhvwmQ5jhDwNCxJii+/eVz71+lVeZKtlSV/z\nHZnlxQg/yXQYI+RRWJAQWyjzVffMNYVOrLPJd5rvxPPng0gEYjE89xzU1Dh4L5nipcgp9l5zrfhP\npnMWo6v5YoqCQYNg3Dh44w1oaXHwXowR6hEWJMQWOrU+PR24Att8Fwy8vywQ6/VQVQVbtsDcudDe\n7uDt6emgU+vpYrX3WgwAfpbpnMVIFMwr+KO4sxPOnoX162HBAsdvxxgh17AgIVYgZ3KysgDAQb4b\nGsOr2yt+/30AgBMnoLDQwR4oCgCALvFqsvOrTNdjjK7vFn/6KQDATz/Br7862APGCLmGBQmxAl2s\nJpmOIPlOa7qV76CV9+Qo8eefg1wOsbEO9kBRIJMBXdzzRGR38bdM5zBGN83dzt3NjekaT2pudrAH\njBFyDQsSYgW6RJ2e3m0JV2CSF9cdPtst393JEc+ZAykpjncik4GussGTzbzFDzOdwxiFTa4rVnU7\nbjj1sTg8HO680/FOMEbIBSxIiHnkjpzkfI61MxdNqzZ2q0mRIbzaPWKDwfF+EhI81EBbfpjpnMUo\nRmIyj+gWo6ExvPOfigMDHe8HY4RcwIKEmCeSCHVqvf3IUEsLXKsxrdpYd63hVr6Li+Bp9zqdZ+yF\nu037Z6ZzFqNjx+DxpaZVG7t9Twrq5NXtwRihPsOChFiBSpWQ23Fai4kBAKg3mMYs6DZ+HsJzfO1L\nXh5IM52cznMTf850DmO0bRu0tIC22TT1qW4x4rZjjFCfYUFCrCBdnELTtgvJCSIuFwJDHVz7Yp/v\nFArw6PMO/DzTOYxRUREAQEYGxEgwRuh2YUFCrEClSmgabA7AH3sMAECvh/feA67ApKiyneNgne/k\n8q79eKiFmOkcxmjkSACAo0ehuhq4ApM2sq7kHMYI9RMWJMQK5KluNkMUUinMnQsAsHo1DB4MCxab\nVm2sO3HZcb4j7xVJhJ5oHmY6cBKj5cuBy4WaGkhOhnvugQlTTS/902lNwhgh17AgIbaQLk6Ry8Hm\npNDXX8NTT4FQCDodAEDiGNPIuQ7OC5Ejd9maVE80DDOdhX2M5s2DL76AkSNBp4ODB6GxEaLEJtGd\nGCPUH/gIc8QWVKpEp07NyCg5ePDW9OLgYPj0U+jogAsXID4ewsMBwGQ21DWdijJpg8k2XcfgsRxP\nPJkUM501hzHKzITHHoMbN+DGDRCLITISMEaof/AbEmILkUQozRxLPZqakWH7PSkgAJKTSTUCcHLf\nmoNrw9zeJMx0NpzFiMOBIUNgwgRSjQAwRqhfsCAhFnFRk2zY5zuh0bD09B43NgYznUMYI+Q5WJCQ\nUwqFAgC8fG9m63y3bJntnC6CjEb8ZYNp3mvd5t31mO/oYrUyX9WbHvlKpiN9wRixGemLwuH/JtQd\nB5/0jqzRNJ2Xl6dQKMjnJ1IUq9VViyRCkTicShNTqRKPXkRioVPr6RK1coeKLlFTVNdNOdPTobAQ\nFIquDCiSCKWZKQ+tks4qPx7fWGt5r54vyBv/gPWulPln6eJKcpPp3vSI5ZmuHz3yUDMwRs447BFF\nURRFyWSy9PR0mUzGdBvZCAsS6rJu3brs7GwAiBTFRorikqjJadJ5kaI4rU5zkT55gT5eotwN/00x\nnhiadohkPXLUTH47SVLWDRAaDQ7znSKnWJFT0tcesTnT9a9HnoYxstbLHlEUlZWVtXbtWqbbyy5Y\nkBAAwLJly3YV7EuVzh9FTU6iJjvcRqvTAECxcneJcte4FZTX8p2FTq13dgmLfb67VmOa9XpHX3vE\n5kxXsHqv+sd6jNFAitHLq5/HmmQNCxLqqkZZC7OdfYRsaHWaHPnzjOQ7F+zzXXNb9B7Vpt68l/Qo\nbAxPJBGyOdNhjAZejLAmWQsgZ2mQ3+prNQIAQXCYdLSsYPtnLa2N3hmu6A0jL7AqLDra0CBs63o6\nRVCAISHy4OXaeT2+l/RIsf+7BNkd7OmRRV8zHWCMvK7fMcr9MLvF2IRDSgQWJFYzGyo5gR65zwrR\nj2pEsDbffb75UnRLw9CYriu++5rv+tEjodFg5Dl5+I879CPTERgjC5bHCGuSBRYk1jEbKk0N59rU\nX7Vc2Nhavj142HIP/SK5XP7+xs39+BQRluzAj+YOSYlxe/P6QZmvKnr3arz4N8Pj2kKCus4L9SPf\n9dgjodEQ1aIbU1cxTVM2oebSr3eMdE8H7CjzVWe20hgji4Eao9wPs2OGREqlUrc3z7fgGBIrmA2V\nHYZKk/ZIe12pSXvEelVY2heBUdM98UsTExNl0qxUaQ9Z4Fr1+S9++DsAPPHgG5IhSTZrdyk2F135\ncnXpc55oYV/lTvt44eTXUqXzBEE1dyZ8EB2qsqwiYxXK84rDJ3deunaKHyQYOyJt4b0vhYVE2OzE\nWY+ERkNYW3O8vlbcWGc9CgIA34yeWRUW7dEe2a/6bNfbVTWX56U/mzLCdpTIZpUPxejJ/yex6VR7\nu3F34cenyg7UN1yPCB8yfcKD96U9vbdI7hMxAidhsv87VBz76iz9Y3l5uSda6EPwXnaMcVGErJm0\nRzxRkORyOU3TSQsnud7M2Nay5eu3btRVAECrscl+g1HU5N2KzXSxmvGTQsp8lU6tT3p0EgAY2mKO\nVfzGOt+FBNW21i/+cEd5Z6eZw+G0tDYVnSxQX7/wPyu2BwR0+xRY98hFgrMWr6/1RLKz7pEN1eXi\nopMFnZ3mxmZdj6t8JUZHzl47fPKYuROsO7X5P2+cvvgLeX1De23ngU0VmrJ7py9hf4zASZhKT+/Z\n/u1am7/Dh2e9vFuxWaFQ+PmJOyxIXtXLItT9LVWea0+kKM71Bvl7c0g1cr6HWACgS5hPdoSlRzb5\nrsXY+act5Z2d5iljZz85782jp/d88cPfKzRlFytOjhk21XoPkhjeXWODp2rK7j5f7iLBWbMM0XuC\nTYy+L9xyRX363JUjnZ1mmy2drWJ/jPK+P3z0vPHgqRZz9/M1V9SnSTWaMfnhtIkP7S/5/IRqv/K8\n4qH0uWyOETiPRVt76xc/rLf/O9Q3aQGgsLAQCxLyrH4UIa5AzJc8youc7qGTdQCQl5fX48m6k+d+\nLjpZEMwPaTU2O9uGXPpHF1fCGnc3sY+UO1Q2PbKuSbuKm/UGMy+As21NXUlFWPrUxwSDhB0dpvDQ\nKAAQBNWEBNVGh6qiw1TRoaqn7xoC0AiNjS5+nZ4vKItMqBJGe+hEkMMeAYDiWL6+qd7h9s5WsT9G\nm7/fr20wWhZOiJcDzAWA8sqzABDIC1o853fhIQ0xwvQTqv0AkDrsn6/OZW+MwHksTp470NLaxOUG\nPPHgG4LgW3+HkiGjkqjJCoXCz6eAY0HyrPa6I43FS/r6LrOhsuXCRoCNnmgS8fVvyX8XOdugqs70\nP+9ouFzO0gVrP8r/vYtdpUnnF5z4h5vb13d0iVq2cIXNQktNOn7xEACMGxb08XdXvvllhr5l0Ohh\nd2bev0YYGhEdqkof2ecsIDQapmnKQFPmntY78vIbAQBnbWLU8CRH3zy4scW88esGAJia8G7ykKZz\n1zMB4MH051paG1uMzfuK8mx2xfIYzUtfERH8s7njOukUn6dPH7m28NK6qePuHzNsaiAvKCqsbFj4\n2h9+bgCAQXyOdHhQb34dUzECuzDFh5f+t8SeAYChsaP2FslPqPa3tbda/g7TpPMVSrnnmuoTsCB5\nkNlQ2Y9qxAbmTnjxn3W6JvPyeXcNl4x3vXESNUlXoGd2iILctCaJcnAqn9Skq9XHARpPXTKeukSO\nxFuOndl35dqvG17b1I9qxKAVc8MAoFbXQTKdNdmdjwGAvqneviCxPEZpk5YLguYNE/1j49f7AcBo\nEhZeWgcAwtBIYWikIKjm8qXfP7616wvHptXRQYEcLza8P5yFqV5/AwDoqnN01TmyhPwdrlv1VRI1\nSV6Q7efDSHi3b+RA7lcNh8+2Th3Nf2Zuz9f5R4riIkWx5CaSzHI2JGZoi6ltjCWv77976Tuv/7Tw\n3pcAoL7h+qETP3mvfczxiRgpK7uucPi1Kstm7R2DA8YkBAVwAQDWfXKzstbkwSZ6Unt7143Pbf4O\ni5W7SIwKbR4R72ewIHkQVyAOmbiB6Vb0x38UTQCgrjW98PdP3/v8FbLwk51/IfO/7Wl11SJJuMNV\n3iESCwHgIn3C2QaCQUMAgB8kWJDxYljI4Dl3Z5EJ31errh+r+I3X2skg9sfI0BZls0Snr7lRV3Hj\nZqBk6O+L3o3bsz4WAMqr278+5HRQk+XIX53932Hl9UsAoNVVJyQkMNxERuEpO8/iSxYFRk5v1x4x\n1ZWSqQ29eRdXIA6MnM6LmsaXOB3juR00TScmJs6TrZwvW+lwg0bjYwDl1dqObBuDFgAAIABJREFU\nau11gOtkYU29OrIu1n5jcvdiaWaKJ5raSyTVXqBPOLs4kRyYczgcLpcLABwOl8cLBAAeL7CiPqO2\nKYXMaCBTG3rzG/V8QVVYdFVYdFmURzKITq3Pnfaxsxjpm+oBZgPA0YpXuKK5Pe7NJ2JkT16QXXb1\n6JSxs59b9HZtU0qk4Exw0FutbaYzV9t683ZmYwRWYapqmEZOlbr4OyQxysrK8kRTfQUWJI/jCsR8\nwSJSWsyGyt4UJ7Oh0mj4yqj+qvnU63zJIrd/zSIPZXFxrProfa8aWvXkdZOhIX9vDgDMnbl8zDAH\ns/4u0MeZzXQAIJIIqVSJix5NTpm155dtrcbmA6U77p2+5PjZH2823AAAyR1JAGBoi6moj6mozwAA\nQVDN5cvrR45Rj58c6WJKsdBoEBorxtRVzCo/XhaVsD9xipd71Cc+ESN7Q6ITy64e/fW84lr1efEd\nSbt/qW1tMwFAwKBF7+6kfTFGLv4OL9DH/bwaARYkL+tHceJFTfNES2Qy2Xu5HzlbOy7pbsvrhsY6\nUpCSh08fmeBgRPoifWJG1gRPNLJPqDTxma20s7WSIUljR9519tLh/L05uxWbDa2NABAdIZ5q9Zg4\nwtAW8/43F2dkT7g8OkVoNMQ31sY31pIXznbuoVnFrnvUJz4RI3v3TltSdLKgvd349uan+UGDWo0G\nAAgViKaMzdyw/VlfjJGLv8PvD21Z+erbnmikD8ExJMZwBWLy7Ud41xeiWb+ETNzAlyziRdp+BQm0\nW+IWS5cu1eqqe3Nwx+G4+iMpUe7W6qqpVOavuJRmjnXdoxcWb0iflBYqCCBZIImatPrpDwJ5thOI\nrXuk5wvIkfU3o2fmjX9gf+KUsqgE+9TmoWTnokcc57PM7Ff5SowsLZ8QLycnTqMjxK88+V5czPDO\nzk5SjUYmTFz99Afny4+xP0bgJEwvLN4wbfzcYH6I9d/h8bM/anXV/jy/jsBvSKzg4psTVyD2xG+k\nKAp6d0JfGBrxUfZx19s4eyabN5E2uOhRnOjSN2s1HWbJGTqsoT2z1mD73ch+b9b0fIGen0BGI2yO\nyvV8gZs64aANDnsUFuI0KM5WsT9GYSER37yTlz5yLYC5ue39c9WZFfUZSdSktS/taG7R39TfiAyP\nHRQcCgCVNy4B62METmIRyAta/sj/ms3mG1paFBZDenSp4hT891Ppz/Bu36zDCRTywpODYmfzhy4a\nNGq1534RRVHb8jZJR8sEwWH92wN5wtjC3DksuZO0SBJ++MtDznqUHJsvGkRzOTBkcFtIcBsZMbLR\nyx4ZeYF1AtHVwXFlUQlH45Pd1gE7rnvUG74YIwAICjAE8QyWGAUF8oWhkeTr7MCIEYfDCQsZbN2j\n7du3492+8ZSd/8rKynp59fM58ufJM5X7SqvT/CH3IWlmCuOj5RbSzJRxKyhnPbKePicIqrHfwOd6\n1COf65HfxigrKwtnNAAWJD+3du3a/tUkS15g28OkZWvSOsKa7XsUHaqyPH0HAOynd7O5R/3Ld2zu\nEcaIsFSj7du3e6htvgULkr/rR01ibV4AgILVe6lUiX12sD/ctk52fe2RTq13S2t7qR/5DmPkKzHC\namQNx5AQyGSyFmNT7ofZhtYmjstnUmh1mp+PfLHpy9+xNtMBwMLcOVSapKW1sWD7Z5YejYz5ngxO\n3MKBivqM3vdIp9Yf2XJCma/6cvnOI1tOKvNV11W1rXpja4PRC3dAsO+Rsy0xRj4UI6xGNvCJsajL\nunXryNFJpCg2iZo8ippCXmh1mmLlbq1OQ64kF0mE0swU2Zqe73HnZZZMZ1miyClW5JQAQKQodv+G\ngKEx3eaUNrdFD33sGPSiR5b9AABFQVYWpKcDTUNeHigU0Js9uIt1jzBG1nwuRhRFZWVl+fnDJuxh\nQULd0DStUCjy8vIU5HP8X+TzTKVKWPKQNxv2mc5Cp9aHnLv8WoSDh0P/+dfI5uThrnskfzQfqtUk\nwdlfJULTQNNQWAhyOYimeukbiU6tp0vUyh0qm5ulYozYHyNSh9LT0/GSI4ewICHHaJqmaTojIyPr\nq0x2JjgLF5mOIDePsV9eGjfG9YRgkunKHeRJWzQNGRney3eETq3XqRvki/IxRuyP0cGDB7EI9Qgn\nNSDHyP3uAMDXMx0AOLujjLixzsW7ep/pAICi4OBB0B1VkfZ4h0giJNHBGPUGszHCatQbWJCQD+tN\npgPnyS6szelTDPqU6QhG8h37YYxQ72FBQr6q95lOaDQ4XOXsdpyKnGK6xDbTdXbCxo1w550QEgKJ\nifDii1Bjd+GmJd+Rp6MiL8fI4uRJSEuDtDRQKm1XYYzYDAsS8km9zHQA4CzTEfF6B8lOma+yv2r+\n1Vdh9Wo4fhwMBqBp+PBDuPdeaGmx3YxM8VLkFPfYsAHP+zEimpvhiSegpARKSkDv6GIkjBFrYUFC\nvqf3mQ6cnwsiHA5R6NT69PRuS8rLYdMmAID774e9e+HxxwEAzp6F//zHwT7T00Gn1tPFzD8vnEHe\nj5HF6tVw4UIPvxFjxE5YkJCP6VOmg56Snf0QBTmTY3P0ffgwmEzA5cL27XD//bBpEwQGAgCcPu1g\nn+SWzTbzff0KIzEivvoKtmwBYU93NscYsRMWJORL+pHpXJ8Osh+ioIvV9pnOaIThw2HaNIiNBQDo\n6ACzGQAgxtH9pikKZDKgiyt72cgBhqkYAYBaDStXQkAAbNvWwy/18xixFhYk5DP6mumgp8EJwmaI\ngi5R258LWrECLl+G4mIAAIMBnnkGOjqAy4X773e8T5kMdJUNvW/ngMFgjMxmePppuHkT/vQnSOvF\nDRn8NkZshgUJ+YZ+ZDrofi7I2SParIcoyB05XTwm7fRpmDoVfvgBAOB//xcmOHkseEJCn5o5QDAb\no7/9DQoLIS0N/vjHXv1e/4wRy2FBQj6gf5kOuic764dYk4eKEtZDFCKJUKfWFxY63tvWrTBtGqhU\nEBoK27fDW2+5+tVevts04xiP0eefAwBUVEBqKsyd27VwxQpYtcrpr/a3GLEfFiTEdreT6axPB+mD\nBA5f2wxRUKmS7rfx65KbC88+C62tMHEinDrleETdIi8P2PMEOS9gSYwAoKoKjh27dfnR5ctOZ9z5\nW4x8AhYkxGr9znTQfXDC+nCbsD4Ytx6ikC5OoWnbXV24AK+/DgAQHQ0ffgjt7VBWBmVlUFXl+Fcr\nFGy/nY8bsSRGGzbAp592/cvN7Vr41ltOz+D5VYx8BRYkxF63k+kAQM8XlEUlkKRmndqIyrAosk1p\n3JhGfohlOZUqoWmwOQDfuhVMJgCA2lqYNg2Sk7v+/e53Dn6vXN61n/4127ewJ0YPPghPPdX1LzOz\na+Hs2eDwciW/ipEP4fW8CUJMuM1MBwBVYdEkxwmNBj1fMLXqnPXa81GUw9tIiyRCKlVSWKi2vhmm\ns9M+XEdHdGR4QyTp6VoY38eqGFlzGBdr/hMj34LfkBAb3X6ms+Zw7pazCV0AIF2cIpeD9UmhnTuh\ns9PBPzKQbo0cucvWpLql5WzGthhZu+OOrgDNnOlgrf/EyOdgQUKs495M1w9UqoR6NDUjw2m+c4am\nYdkygFgJC5/W6l4YI+QJWJAQuzCe6aDr0atj+5rvSKajjZKsrzN73tqXYYyQh2BBQizChkxH9DXf\n+U+mwxghz8FJDYgt2JPpCJLvACAjo0Qmg6VLwX4InaaBpqGwEORygNiBn+kwRsijsCAhVmBbpiNI\nvhNJwhU7VPIMNUV13ZQzPR0KC0Gh6Jp5LJIIpZkpA35MAmOEPA0LEmIeOzMdIZIIpZIUaWaKTq2n\nS9R0sTo7W0WWi8ThsjViP8lxGCPkBViQEMPYnOmsWbLewtw5OrXery5hwRgh78BJDYhJvpLpbPhV\npsMYIa/BgoQY46OZzq9gjJA3YUFCzMBMx34YI+RlWJAQAzDTsR/GCHkfFiTkbZjp2A9jhBiBBQl5\nFWY69sMYIaZgQULeg5mO/TBGiEFYkJCXYKZjP4wRYhYWJOQNmOnYD2OEGIcFCXkcZjr2wxghNsCC\nhDwLMx37YYwQS2BBQh6EmY79MEaIPbAgIU/BTMd+GCPEKliQkEdgpmM/jBFiGyxIyP0w07Efxgix\nEBYk5GaY6dgPY4TYCQsScifMdOyHMUKshQUJuQ1mOvbDGCE2w4KE3AMzHfthjBDLYUFCboCZjv0w\nRoj9eEw3APk8lmc6nVqvzD9LF1caEhqnLRGRhcp8lbzoLJUmplIlVJqE2RZ6AcYI+QQsSOi2sDnT\nKXKKFTklABApio0UxY1KTwK4QlbdETAytmWQImc3QIlIIpRmpsjWpDHaWA/CGCFfgQUJ9ZnZZFbm\nqypPVF/cfyVwUOBdv5na0W4OCGTX6d+C1XvVP9bPk60cRU1OoiYDQNKQfEuyS6ImZy3MnC9bCQDF\nyt0lW3cBwIDJd22G9l/eLSWvL/50FQCS7hv28/8rEk+KHTV7OKNN68afYwQAtbW1mzZtOnnypE6n\nS05Onjlz5uLFizkcDtPtYhIWJNQ3zXWGHc9+d+1olWXJ7v/5qXTLiWXfPi6IGMRgw6yRTJe1MJuk\nOWciRXEAMF+2Mk06L2fr8zBQ8l1NWZ2lIBE3ymoBYNKScewpSH4eo6Kiojlz5jQ3N5MfCwsLN23a\ntG3btl27dvH5fGbbxiB2HdUi9vvpr4dINeJwOeFxYWRh7aX6na/tY7Rdt/Qy01mLFMWtyfrozFZa\nkVPs0bZ5Byk/ABAQFCAYPMjyjx8WxGzDLPw8Rh0dHY888gipRgsWLHj99deTkpIA4Kefflq/fj3T\nrWMSfkNCfaCvbvz1P+cAICAw4MWfn4kaEXHhpytfLC0AgMsHy01GE4/P8F9UPzIdQfLdwDgGv1FW\nBwAcDrx5YRXjEbGHMTp9+nRtbS0ALFu2bNu2bQDwxhtvDB06tLm5edeuXX/+85+ZbiBj8BsS6oPa\nS/Wd5k4AmPh4StSICAAYdd/wCEoEAB3t5nq6gdnmKfNV/ct0hOUYXJmvcnvbvOnc7gsAEDk8gsPl\n3Cir1Zy+0dHewXSjumCMAMByUi4oqOs7a2BgIBk9ioyMZKxZLMC6oyfEZoc/ODpIFBweL0yYLiZL\n2gzthvoWAOBwOYOHChltHShyihfKXnOW6QpPXXwtt5q8NrR91tK+FwBWPZEbMuhWsyNFcanS+Yqc\nL6WZKV5osCcUrN7b0mAEgHZDe+7UjxtvNANAQGDA3avulK1J43AZHjN3HaMNn++rrq1+fbHo3kld\n45Gd0HngyJelp/dU114NC4lIGZH6UMbzvh6j5OTkESNGXL58+dNPP01KSho9erRcLm9qagKA+fPn\nM906JmFBQg7QNJ2XlyeXywEgd9rHZMZtweq9wtiwZ3Y8ZtmsvdX0zaofWvVGAEi8a2jgoECmGgwA\nynyVTq1PenSS0w0uqo+dN/73Jw2ABgA6Okw2m42iJu9WbKaL1Sy/9oVcu0O+KFjHqL2lvaOtAwAa\nNI2WjTvaOwr/eYTD5crWpDLW4p5ipLpcvKvodGdnp7bh1ve5HXveOVi6g7zW6jSHjn99+Zrysdm/\n9bkYJSYmZmVlrV27lqw6dOhQampqRUXFmjVrLNv/9re/feGFF5hpKztgQUK2aJpOTEy0/KhT6xU5\nJYqcEmlmivW1LNVnawpe3UOGK0IiB83/+30MtNUOmZTl0FVNLQBk3R82a/Ig9c271DfvBgDBoDC7\nPcQCAF3C6mSnU+tzp31s/aMlRtLFY1W7LgJA7NiYhRvnhA0JLd16qvAfJQBQvOlY6vOT+aEMT22w\nj9H3hVuuqE+fu3Kks7PTenndTU3hsa8AIGVE6r3TnyhR7jp29kdNzRVN7RXwtRjRNJ2dnS2Xyw8e\nPCgUCufOnVtRUUFWhYWFNTY2AsD27dvHjRu3bNkyZlrMAjiGhGw5+zzQJWryotPcefhfx7Y8+Dmp\nRkNSYp79/snBCeHea6Ijyh2qVOk8Fxtc1dQBwEN3Ce6bMujO5MQJo9MnjE7nBdh+q4sUxSVRk+ni\nSg+29baRa13t0SXq4TMT3rz48psXX37uhyfvGBMtGDwo43dp8dIhANBmaK+7VO/dlnbjLEaKY/mq\ny8WdnWab5VfUSrO5g8PhLl2wNmVE6hPz3gwI4AGArrHWR2NE03RGRsYXX3yhVCoBIC0tTaVS6fX6\nffv2BQcH63S6VatWtba2er2xbIEFCXUjl8sVCoXDVTq1ni5Wd7Sbv3rp+5/+71BHu5kbwJGtSX3u\nhycHD2W4GgEAXaIeRU1xtrax+aau0QAAHxTohy6+NvvVf/7t42cqr190uHGadL6ukuEJGi4o81WW\ngwMbOrW+bO+lhkp9Q6Xe3HErv5OJJwDQXGfwRhOdcBajB9OfW3jvS/ffvdRmucnUHh0hThSPDQ+L\nAoBOs5kUrbCQCN+NEU3T5GQ4ALz55pvJyckAMHv27McffxwADAbDzz//7K1msg6eskN9oMgpadUb\nr6tqAIDH5417eDQA/PLuEabbBTq1HgCSKKcDSJqaruv/fz7Zwg/kmM2ddNW5d+Qr//rKzhCBbTVN\noibpCvQFq/eKJAxP03CIdNaZH9480Kw1AMCwGUOHTo0HAHNH58Wfr5K1V4sqNKeve6GR9lzESHbn\nYwCgb6rfV5RnvfyuSQvumrSAvG5rb9327Z/NZjOHw00ZPn1QcKjvxojM+QaA0NBQy0IyqQEATCbb\ncU3/gQUJdWM5r+1QoCDQctwXGi2ou1Jfd+XWKaB46ZBgIZMXmbsYQGpu0cdGD+MF8J6Y92aieGzR\niW8/2/V2S2tT0ckC+wPzSFFcpChWp9b7YrITDQ0nBenaUU1nJ4TGhFSeqDY2tgFASJSA2QCByxhZ\nHK14hXtqrvWSyhuXtn79lqbmKgA8lPGCeEgSAPhujMaNG0c+aC+99FJeXt6oUaP+/e9/f/311wDA\n4XDS0nz4EqvbhAUJdZOenu5ibUjkrZsD6Sr1uspunzrZa6nD0ykPNcw1ulitzFddpE84m088Kfme\nScn3WH6cMfmRHw9/WlOvrqq57HB7ra56xuI57JxYTKeq5Yscnw4CgIzX0w6sL9L8esNkNJUXXbMs\n5wUFPLNj0R1jor3SRgd6jJEzh0/u/OKH9e2mNn6Q4PG5r6dJuyZG+26MVqxYUVZWduXKlbKysqlT\np1qv+stf/hIdzViMGIdjSKgbmUwmk8kcrqJSJWQ+MQuJJOEAcIE+4WyDfUV5n+36655ftluWmDvN\nABBgN6kBAEqUuwGAnZkOAKg0CZXqeHYZlSoZIaOe/OzRqcukvKAAy/KhU+NX7nuawWoEvYiRQz8f\n+fcn3/2l3dQmiR31xxc+t1Qj342RTCZbuHDhgQMHli1bZrkwFgBGjBixbdu2t956y1ttZCP8hoRs\nbd++PSMjg6Zp64UiiVC2JpVKkzz6rwcZapcrIomQSpVcdJ7sam9W/nLiWw6HQ8WnjBk29ZcT39Td\nrAKAxHgHGe0CfZy1mY4QSYRQ4mAhucwoJHLQ3L/ee3+27Oa1hladMWpkBONn6qAXMbJ3o67iqx83\nAkBYyOCn5v2ho8NUXVsOAIP4IT4aI4qiyKVIQ4cO3bZt27/+9S+apnU6HUVRQ4YMYaCVLIMFCdmi\nKOrgwYPkwliapn3lUTRUmvjMVtrZ2rsmLihR7jZ1tOd+8lIwX9BqNABAbPSwVKmDC+Mv0idmZE3w\nXFNvE5lPvLr0OXLRJRlHsY9RQGBA1PAIhtromOsY2Ss6tdNs7gCAxuabf/v41lDfnWNnX60841sx\noijK+sJYIjg4ePTo0Qy1kY3wlB1ygBzHlZeXA8Dq0uesM93UqnNTq849cv6Q0MjkBGJ70syxWl21\nswPwRPHYFx/PiY0eBgCtRgOHw5kwauarT78fyLO9SrREuVurq3Z2voVxlqftiSRC2Zq01aXPwUCJ\nkf2TgG7U0Q631OqqfS5G5eXlNtUI2cNvSKhvpmnKyIuHLxSWRSYcjU9mtj0WZLbVBedj5mNHpo0d\nmdbYXN/QWBc1WBzMF/S4N7bp5bNffTRGYSERH2Uft17y0pJ/ONxPiXL31cozPh0j5BB+Q0J9YH3E\nLTQapmnKlp7eM7XqHINNsrYwd06JcpdWp3GxTVhIhHhIkrNqpNVp5AXZ7Mwmvcx0AyBGrg2AGCFn\nsCChvtF3T+WWlMeGs0PSzJRxK6gc+fP9y3daneYPuQ9JM1NYOFrep0yHMWIEVqPbhwUJ9YGeLzga\n5+D8j9BoePhCIRsOw2Vr0vqX7yyZjoUJpa/VCGPkfViN3AILEuqbsqiEqjAHl7Ow5+xQP/LdAMt0\nGCMvw2rkLliQUJ/tT3R6D1OWnB2y5Ltdis2uL3zR6jS7FJsHXqbDGHkNViM3wll2qM/0fEFp3BjL\nVC575OwQs/O7yDTo3TmbdwNEimKTqMmjqCnkhVanKVbu1uo05Gp/cj0pCy+0up1MhzHyDqxG7sWx\neSIWQtY4HE62Zo39cpLOejzE1vMFjE871qn1dIlaucP2cQDkYlIqVcLOh7z1PtNlx+VgjBjRpxhh\npu0NLEjIFWcFCQDG1FXMKj/ucJUNPV/w7ah0vcvrfrxAp9br1F1P0GFngrPo03G3s4IEGCNP6muM\nMNP2Bp6yQ/1UFpUwpq4ivrG2xy3ZcHYIAEQSITsvpbThxrNAGCMPwTN1HoKTGlD/uRg5t8Ge+V0s\n5/ZMhzFyO6xGnoMFCfUfGTl3sYHN5GOS8sbUuXoGoD/zRKbDGLkXViOPwoKEbsv5KMrFwENYW/P+\nxCnWG1SFRZdFJXilaT7Gc5kOY+QuWI08DQsSui3O7gtACI2GqZpz345Ktxykl8a7Olr3Wx7NdBgj\nt8Bq5AVYkNDtsr8vgPWPZLT8aHxy3vgHSuPGOLyDgJ/zQqbDGN0mrEbegQUJuYH9yLlNvlt6eo+e\nL2DPcxDYw2uZDmPUb1iNvAYLEnIDm5HzsLZmm0Nyku+YaBqreTPTYYz6B6uRN2FBQu5hPXIuNBrG\n1FXsT5yC+c4F72c6jFFfYTXyMixIyD1sRs7jG2vjG2sx3znDSKbDGPUJViPvw4KE3MbmFNBUzTk9\nX4D5zh6DmQ5j1EtYjRiBBQm5k2XknNyyk7zAfGeN8UyHMeoR4zHyW1iQkDuRkfPSuDF54x+wzNfC\nfGfBhkyHMXKNDTHyW1iQkJsdjU+2nzqM+Q7YlOkwRs6wJ0b+CQsS8hI/z3c+kekwRsD6GA1sWJCQ\n9/htvvOhTIcxYrohfg0LEvIqP8x3PpfpMEaIKViQkLf5Vb7z0UyHMUKMwIKEGNC/fKdT6z3cLjfz\n6UyHMULeh48wR8wg+W5W+XHLA7ZJvssb/4D1Zjq1Xpl/VqfWK/NVACCSCKlUCZUmEYmFVJqEgXb3\n2gDIdBgj5GWczs5OptuA2IvD4WRr1nhu/0KjwTrfAYCeLyD5TpFTrMgpIQspCrKyID0daBry8kCh\nAAAQSYTSzBTZmjTPNa/fvJnpsuNyMEb94OUYYabtDSxIyBVPFyRwku8y1jVCtZokOJnM9i00DTQN\nhYUgl4NoagrbjnC9fNzt6YIEGKPbhgWpl7AgIVe8UJDAUb67VmOSPlvZ4xtpGjIy2JXvvH8WyAsF\nCTBGtwcLUi/hpAbEPPvx86ExPN1+cY9vpCg4eBB0R1UkxTBuAI9JYIyQF2BBQqxA8p1KF2RZYjb4\nWL4b8JkOY4Q8DQsSYovv3lc+9fpVXmSrZUlf8x2Z5cUIP8l0GCPkUViQEFso81X3zDWFTqyzyXdl\nn4rT0sD6X3297XvJFC9FTrFXW/xf/pPpeh+jvXbfhTBGqEdYkBBb6NT69HTgCmzz3R3hvH8tF5eU\ngOVfe7uDt6eng06tp4vV3msxAPhZput9jGprHbwdY4Rcw4KEWIGcycnKAgAH+W5oDK9qp7igAMi/\nwYMd7IGiAADoEq8mO7/KdD3G6MoOMZn/zeHAyJEO9oAxQq5hQUKsQBerSaYj7PNdMPBmCsQLFsCC\nBRAUZL8DoCiQyYAu7nkisrv4W6brMUaiYN6/lokBIDsbpk93sAeMEXINCxJiBbpEnZ7ebQlXYDIm\n1B0+eyvfdbbwyj4R//qr053IZKCrbPBYG7vxw0znMEa25+5EvLJPxH/6k9OdYIyQC1iQEPPIHTnJ\n+RxrZy+ZVm3sVpPuEPE4J8RareP9JCR4qIG2/DDTOYuRfU2KCec1/Ox03h3GCLmABQkxTyQR6tT6\nwkLb5fX1EBpl+tdPdYbAW/lOHMVrUjjNd16427R/ZjpnMYL/1qQKXW/ngmOMkDNYkBArUKkScjtO\na48+CioV7DloikvvdgweFug43+XlgTQzxZPN9OtM5zBGBFdgem5Dt++yzmoSxgi5gAUJsYJ0cQpN\n2y5cvx5WroS//c3BeSGH+U6hAI8+78DPM53DGBFXr8KJs6ZVG+uaeRgj1H9YkBArUKkSmgabA/Ar\nV+Djj+Gtt2D/fuAKTN+dd3UMLpd37cdDLcRM5zBGxC+/AADcaDDF3O3quAFjhFzDgoRYgTzVzWaI\nYsUKCAqCzk647z4QCuHpZ02rNtaduuo435H3iiRCTzQPMx04iRFRVAQAIBYDX+jquyzGCLmGBQmx\nhXRxilwO1ieFpk2DggJITgYAaGwEDgek002jH3KQ78iRu2xNqicahpnOwj5GREkJAEB8PICjeXcY\nI9RLWJAQW1CpEurR1IyMbvnugQdApYIbN0CphIYG2LkT4oc5OQaPlXjiyaSY6aw5jBEAnD0LnZ1g\n+fLkdMwPY4RcwoKE2EIkEUozxzrMdzExMGEChIV1/ejwvjUH14aBu2Gms+EiRjYwRqgfsCAhFrmd\nfCc0Gpae3uPGxmCmcwhjhDwHCxJySqFQAICX781sne+WLXM8p4uMRvxlg2nea93m3fWY7+hitTJf\n1Zse+UqmI33BGLEZ6YvC2TVcyAoHn/SOrNE0nZeXp1AoyOcnUhSr1VUXd3ADAAAgAElEQVSLJEKR\nOJxKE1OpEo9eRGKhU+vpErVyh4ouUVNU100509OhsBAUiq4MKJIIpZkpD62Szio/Ht9462kHer4g\nb/wD1rtS5p+liyvJTaZ70yOWZ7p+9MhDzcAYOeOwRxRFURQlk8nS09Nl5L7oqDssSKjLunXrsrOz\nASBSFBspikuiJqdJ50WK4rQ6zUX65AX6eIlyN/w3xXhiaNohkvXIUTP57SRJWTdAaDQ4zHeKnGJF\nTklfe8TmTNe/HnkaxshaL3tEUVRWVtbatWuZbi+7YEFCAADLli3bVbAvVTp/FDU5iZrscButTgMA\nxcrdJcpd41ZQXst3Fjq13tklLPb57lqNadbrHX3tEZszXcHqveof6zFGAylGL69+HmuSNSxIqKsa\nZS3MdvYRsqHVaXLkzzOS71ywz3fNbdF7VJt6817So7AxPJFEyOZMhzEaeDHCmmQtgJylQX6rr9UI\nAATBYdLRsoLtn7W0NnpnuKI3jLzAqrDoaEODsM1AlgQFGBIiD16undfje0mPFPu/S5DdwZ4eWfQ1\n0wHGyOv6HaPcD7NbjE04pERgQWI1s6GSE+iR+6wQ/ahGBGvz3eebL0W3NAyN4ZElfc13/eiR0Ggw\n8gL72eJe6EemIzBGFiyPEdYkCyxIrGM2VJoazrWpv2q5sLG1fHvwsOUe+kVyufz9jZv78SkiLNmB\nH80dkhLj9ub1gzJfVfTu1Xjxb4bHtYUEdZ0X6ke+67FHQqMhqkU3pq5imqZsQs2lX+8Y6Z4O2FHm\nq85spTFGFgM1RrkfZscMiZRKpW5vnm/BMSRWMBsqOwyVJu2R9rpSk/aI9aqwtC8Co6Z74pcmJibK\npFmp0h6ywGe73q6quTwv/dmUEQ5GI3YpNhdd+XJ16XOeaGFf5U77eOHk11Kl8wRBNXcmfBAdqrKs\nam6L/k6Zu7vw41NlB+obrkeED5k+4cH70p4O5AXZ7MRZj4RGQ1hbc7y+VtxYZz0KAgDfjJ5ZFRbt\n0R45XNtk0P1waKvyvELfpI2OkIxLunt++srAQL7NZj4Uoyf/n8Tmj81sNn9/aMuxs/vqG67HREhm\nTHo4Y9piH4qR/cfnVNmBH4s/s9ls5NCJF6oKy8vLPdFCH8JjugH+y0URsmbSHvFEQZLL5TRNJy2c\n5Hoz1eXiopMFnZ3mxmadww1GUZN3KzbTxWrGTwop81U6tT7p0UkAYGiLOVbxG+t8FxJUW7D3wf0n\nunpxQ3tt54FNFZqyFx9/x2Y/1j1ykeCsxetrPZHsrHtkz2Rqy/30N+rqCwDA4wVpaq5oaq5cqz6/\n+ukPbLb0lRgdOXvt8Mlj5k6w/mPL25l95NcfyOuqG5e/3LOhoakuefh0n4iRw4/PRfrkVfVpmy1n\npz2973CeQqHw8xN3WJC8qpdFqPtbqjzXnkhRnLNV3xduuaI+fe7Kkc5Os8s9xAIAXcJ8siMsPbLJ\nd0fPG0k1mjH54bSJD+0v+fyEar/yvEKrqyZdsJDE8O4aGzxVU3b3+XIXCc6aZYjeE5zF6Ir6tLr6\nAofDeWnJP8YnzSg9/cO2b/5cdqW08vpF8ZCk7ntge4zyvj989Lzx4KkWc/fzNTe010pP7wWAmVMe\nyZi6eOeBTcrzigOlX86QSlkeIxcfH03NFQCYOeWRsSPusiwcEp0IAIWFhViQkGf1owhxBWK+5FFe\n5HQPnawDgLy8PNcn6xTH8vVN9T3uh1z6RxdXwhr3Na5flDtUNj2yrkknLhoBgB/E2faa+sClcQ/M\nWH5CtR8AdI21kaJYQVBNSFBtdKgqOkwVHap6+q4hAI3Q2Oji1+n5grLIhCphtIdOBDnskbXmlgYA\nCOTxU4anAsD4pJlkubG91WZL9sdo8/f7tQ1Gy8IJ8XKAuQBwquxAZ6c5KDB40ezfDg5tfGnRgpX/\npzC2tYgC/vLdX1kdIxcfn6qaKwAwKXnWKGpyR4fJcoo1iZqsUCj8fAo4FiTPaq870li8pK/vMhsq\nWy5sBNjoiSYRX/+W/HeRsw0anuTomweX1yZ9vq/U9a7SpPMLTvzDnY3rF7pELVu4wmahpSY9OuN0\n+oTg4EBOSFBtQujv8o6EAUBQYHBC7OjoUFX6yD5nAaHRME1TBpoy97TekZffCAA46yxGsmEdW7/m\ntLW3HixelUi9dPhkAQAMCg6l4pLtN2Z5jOalr4gI/tnccX3j1w0AwOfp00euLby0rr7hOgDExQwX\nR1wlMfpzRMD1+o7KWlNvfh2DMSIfn8YWM+lRfHgpKbGNzTcbm+sB4Kfiz97/92qz2TQ0dvTT8/8o\nHpKUJp2vUMo911SfgDdX9SCzobIf1YglVswN++1j4Zn3TulxyyRqkk6t9/L9PW2Qm9YkUQ5O5ZOa\nxAkcn5wQNCwu8MPdnXN+pyhW7gKAZQ//r1Cg60c1YoOo8ICPXosCgP/8fPzvW5cfPvUdl8t9YfGG\ngAAHR5ksj1HapOWTpP/MvPdO8qPRJCy8tA4AbjbcAIDw0EGWGEUKAwCgqq7DO23uN/LxefEh22s2\nyPk6AFBdLgYAs9lMV517R76y2dCQRE2iadrP78GKBQndrkhRXKQoltxEkvGWOFxOalJtU0pzW7Sm\n+YW4mOFcLhcAvt3/nlbXqxEIFjp1yfibjXUAEMgLGJkwMZgfYjabP9v119r6SvuNfSJGysquKxx+\nrcoiL8ydZgAwm22HYUwdvjo3uLlFHxs9TDIk6X+e3f7+Hw8/Nf8PANDS2lR0soDEqNDhI+L9BhYk\nD+IKxCETNzDdCm/Q6qpFknAGGyASCwHgIn3C2QaaOth1fGH+0bfvHDt77Us7fr98GwDU1KsPnTx+\nrOI33muo+3y+v6m1rVMQzPn6by/8btnHf37xi2B+SG19ZekZxw93YH+MDG1Rtm8JiwaA5pZ2S4x0\nTWYAiI/y1bGGScn3ZP8m/48v/HuYeBwHODMmPxITIQGAqprLAKDVVSckJDDdRib5alx9BV+yKDBy\nerv2iKmulExt6M27uAJxYOR0XtQ0vsTpGM/toGk6MTFxnmzlfNlK11vqm+oBbKcR2yB3L5Zmprit\nfX1HUu0F+oSzixPlBdllV49OGTv7uUVvA0BC3JigwOC29tZr1y9U1C+tbUohMxrI1Ibe/EY9X1AV\nFl0VFl0W5ZEMolPrc6d97CJG+04uBmgcMfTu6uZlABApipMMSbpUcerKtV/tN/aJGNkbLIwBgOt1\n9NXaGbVNKYO4xzXaPwNAXFRAb97OeIz0TfUAswGgqmEameC4ryiv9mZlpCjugRnLyDbkW2BAQCCJ\nUVZWliea6iuwIHkcVyDmCxaR0mI2VPamOJkNlUbDV0b1V82nXudLFrn9axZ5KIuLY9U+uUAfZzbT\nAYBIIqRSJS56NCQ6sezq0V/PK65VnxffkaQ49p+29lYAiI8ZAQCGtpiK+piK+gwAEATVXL68fuQY\n9fjJkS6mFAuNBqGxYkxdxazy42VRCfsTex5sc2+P7ogcqqm5Ul6lamisCw+Lqq2vvFZ9HgDucJR8\nfSJG9kYPm7pLsdnQ2rj/yL/vnfa4fP/5zk7gcrkm/vp3d37K/hjZq71Z+cuJbzkcDhWfMmbY1F9O\nfFN3swoAEuNTLtDH/bwaARYkL+tHceJFTfNES2Qy2Xu5H7llVxfpEzOyJrhlV7eDShOf2Uo7W3vv\ntCVFJwva241vb36aHzSo1WgAgFCBKHXCgzZbGtpi3v/m4ozsCZdHpwiNhvjG2vjGWvLC2c49NKvY\ndY/SJj6kPF/Y2Fyf/cFjkthRFZpzxraWwEB+6gQHs5B9Ikb2RgyVjkyYeKni1Dc/vfvdgU2mjnYA\nmJxyX1iY9G9b/8T+GNm7a+KCEuVuU0d77icvBfMF5O8wNnpYqnT+3iL5ylff9kQjfQiOITGGKxCT\nbz/Cu74QzfolZOIGvmQRL9L2wqNAuyVusXTpUq2uuseDOw6nh/2UKHdrddVUKvNXXEozx7roUXSE\n+JUn3xs6RNzZ2UmywMiEiauf/mBw+B02W1r3SM8XkCPrb0bPzBv/wP7EKWVRCfapzUPJznWPxifN\neHbRX2MihxpaGy+UH281GobGjn75iY0JcWNc9IhZrntk+WObEC+3nDh9PnP92JFpAGDqaOdwuNMn\nzH1q3lu+EiP7j0+ieOyLj+fERg8DgFajgcPhTBg189Wn3z9+9ketrtrPr4oF/IbEEi6+OXEFYk/8\nRoqioBcn9MNCIj7KPt7j3pw9k82bSBtc9Oiusfw/PMy72TT0SrWw2by41vCAw82s92ZNzxfo+Qlk\nNMLmqFzPF7ipEw7a4KJHU1Lum5wyq7GpXt+kHSy8I0Tgas4C+2MUFhLxzTt56SPXApib294/V51Z\nUZ8RFhLx8pPvthqbtbrq6AhxUGCwzd6ssS1GDj8+Y0emjR2Z1thc39BYFzVYHGzVMPKp9Gd4t2/W\n4QQKeeHJQbGz+UMXDRq12nO/iKKobXmbpKNlguCw/u2BPGFsYe4cltxJWiQJP/zlIWc9So7NFw2i\nBwVx4iLbQoLbyIiRjV72yMgLrBOIrg6OK4tKOBrv4EJUd3HdIwDgAIcfJBCGRlpnamu+GCMACAow\nBPEMlhjxeEHC0AhyiZXPxcghftAgYWgkjxcI/+3R9u3b8W7feMrOf2VlZb28+vkc+fPkmcp9pdVp\n/pD7kDQzhfHRcgtpZsq4FZSzHllPnxME1dhv4HM96pHP9chvY5SVlYUzGgALkp9bu3Zt/2qSJS+w\n7WHSsjVpHWHN9j2KDlVZnr4DAPbTu9nco/7lOzb3CGNEWKrR9u3bPdQ234IFyd/1oyaxNi8AQMHq\nvVSqxD472B9uWye7vvZIp9a7pbW91I98hzHylRhhNbKGY0gIZDJZi7Ep98NsQ2sTx+UzKbQ6zc9H\nvtj05e9Ym+kAYGHuHCpN0tLaWLD9M0uPRsZ8TwYnbuFARX1G73ukU+uPbDmhzFd9uXznkS0nlfmq\n66raVr2xtcHohTsg2PfI2ZYYIx+KEVYjG/jEWNRl3bp15OgkUhSbRE0eRU0hL7Q6TbFyt1anIVeS\niyRCaWaKbI2Dp8cyy5LpLEsUOcWKnBIAiBTF7t8QMDSm25zS5rbooY8dg170yLIfAKAoyMqC9HSg\nacjLA3InTK/9P7HuEcbIms/FiKKorKwsP3/YhD0sSKgbcr/hvLw8m7sOk88zlSphyUPebNhnOgud\nWh9y7vJrEQ4eDv3nXyObk4e77pH80XyoVpMEZ3+VCE0DTUNhIcjlIJrqpW8kOrWeLlErd6hsbpaK\nMWJ/jEgdSk9Px0uOHMKChByjaZqm6YyMjKyvMtmZ4CxcZDqC3DzGfnlp3BjXE4JJpit3kCdt0TRk\nZHgv3xE6tV6nbpAvyscYsT9GBw8exCLUI5zUgBwj97sDAF/PdADg7I4y4sY6F+/qfaYDAIqCgwdB\nd1RF2uMdIomQRAdj1BvMxgirUW9gQUI+rDeZDpwnu7C2/9/e3Uc1deZ5AP8iKkoXCI1vUHK42BW0\napeKpxb2dAltbe34Mu62Mjp9Me6eatdjFzvOtB3tS5zaTru7zMHO2OppXbDrS4v1iKu7cjxjG6wl\n1WpFFF/w7dIgIgjGKAgKyf5xHSaEJPKW3CfJ93P8A3IfksfzPfn9knufe2+Tpz/pUaVTqFLvxMeM\nqPvYkChQdb/SRbc2u93k6XKcptxS2eyx0i1ahIwMFLsraB31Trk7Kvk5o5s38dvfIiUFQ4ciJQWr\nVqGlxfUPmZHI2JAoIHWz0gHwVOkU99ncFLuywgpPZ80XF+Ozz2A2o97DdaWVJV6m3NK7Tizo+T+j\n7Gx88AEqK9HSgspKvPUW5s1z84TMSFhsSBR4ul/p4HlfkMLtIQqrxZaZ6frgu+/i6acxfTq63FDb\nVWYmrBabXKr+/cJV5P+MSkuxaxcAvPQSzGbMmQMARUWoqnLznMxITGxIFGB6VOlwt2LX9RCFsien\n66fvNWtQXHz3bgRAuWSzy3rfkKJKRgcOAMCQIVi9Go88guXL7zx+8aKb52RGYmJDokDSi0rnfXdQ\n10MUcqnF7f66t9/Ge+/h9dfv/qKSBL0ecml1NycZZNTKaN48lJejvBxhYaisxOrVABAZiTR39+4I\n8YyExYZEAaOnlQ53OzihcDlEIZstXffXAVi8GMuX49VXu/W6ej2s1de6NTS4qJjRqFGYOBFjxuCT\nT5CSgoICAPj8c0REuH/OkM1IZGxIFBh6UenQeV+Qp1u0OR+iUK7I2ffbpCUm9vUZApEgGcXFYcIE\nhIcDwBtv4Kef3A8LzYwEx4ZEAaB3lQ6di53zTayVm4oqnA9RaHTRVoutpKSX83Tm56tNq071jC5e\nxOnTqK/H3Lk4dgzffQcAZ89iyxaPLx1qGYmPDYlE15dK57w7yDY40u3PLocopHRd58v49caGDRDn\nDnJ+IEJGBgPGjsUrr9z5NS0NkZEAcOSI+5cOtYwCAhsSCa3XlQ6dD044f9xWOH8Ydz5EkfqL8bLc\ni1frxGQS/XI+/UiQjMaNA4AdO/Djj2hvx8cfo7kZACZOdP/SIZVRoGBDInH1pdIBsEVEnhyWqBQ1\n59KmqI4apow5ED/uesQ9HY9L6TpZRl++JCmH06X0kCh24mSUk4OhQ9HSgsmTERuLnBwAGDYML77o\n5nVDKqMAMvDuQ4jU0MdKB+Bi1HClxkW3NtsiIh++eMJ566lhktvLSGt00VK6rqTE4vZimGFhd39d\n5fCGRhfdizkHFqEyuv9+7N6NJUtw/DiuXweARx/F6tXQuWs6oZNRYOE3JBJR3yudM7drtzwt6AKQ\n+ovxBQVwu+NuxAg4HHA48MIL7v9W+eSuX5beu6kGEAEzyszEsWNoaMDRo7BasW8fHnrIzd+GTkYB\nhw2JhNO/la4XpHSd9Ex6Vpb7nuSFLGPBAiBOJ+DdWvuXyBndey8efBAxHm5ZHjoZBSI2JBKL6pUO\nd269OqGnPUmpdHKrzrAt24eTEwAzIh9hQyKBiFDpFD2td6FT6ZgR+Q4XNZAoxKl0CqXeAcjKMuv1\nmD8fXZc5yDJkGSUlKCgA4oK/0jEj8ik2JBKCaJVOodQ7jS7G9GVFQZZFku5clDMzEyUlMJnurDzW\n6KJTs8cH/TEJZkS+xoZE6hOz0ik0uuhU3fjU7PFWi002W+RSi9FYoTyuSYjRL0sIkRrHjMgP2JBI\nZSJXOmcdVW923jSrxRZSp7AwI/IPLmogNQVKpXMRUpWOGZHfsCGRagK00oUUZkT+xIZE6mClEx8z\nIj9jQyIVsNKJjxmR/7Ehkb+x0omPGZEq2JDIr1jpxMeMSC1sSOQ/rHTiY0akIjYk8hNWOvExI1IX\nGxL5Ayud+JgRqY4NiXyOlU58zIhEwIZEvsVKJz5mRIJgQyIfYqUTHzMicbAhka+w0omPGZFQ2JDI\nJ1jpxMeMSDRsSNT/WOnEx4xIQGxI1M9Y6cTHjEhMbEjUn1jpxMeMSFhsSNRvWOnEx4xIZGxI1D9Y\n6cTHjEhwbEjUD1jpxMeMSHwD1Z4ABTzBK53VYisrPC6XVjcnXp8yT6M8WFZYUbD/uJSRIKXrpAyd\nujP0A2ZEAYENifpE5Epnyi015ZoBaDVxWk18SmYycE7ZNDJ8TNzNoabcXYBZo4tOzR6vX5ah6mR9\niBlRoGBDot5bN22j/bb9KWOm2hNxo2hpsWVP4wz9whQpLVlKA5A8qrCj2CVLaYbZ2TP1CwGUlu0y\nr98JICjrXeHCnQ3nriZPHX2q+OzYaX+r9nQ6CdmMjhw58tVXX7ndFBMT89prr/l5PuJgQ6Je2vzi\n9kvllwGc2Xth9KOJak+nE6XSGWYblTLniVYTD2CmfmFG6ozc9YsQLPWuw/ac3VXm6qaG5ssn69Oe\ne1CohhTKGR07duz99993uyk2NjaUGxIXNVCP2dvsG7K3nt//k9oTca+blc6ZVhO/zLDu2HrZlFvq\n07n5U9HS4kvH6poamtWeiBvMiNziNyTqme05u8u3nXTYHWpPxL1eVDqFUu+C5jN40dJi26Ub9aev\nqD0RN5jR9OnTDx486PzIhg0b1qxZA2DhwoUqTUoIbEjUMxf2W4TtRmWFFb2rdIqOeqfRxaRmj+/3\n6flN0dLi2zfbaivqHA4MT9bWVzaoPaO/YkYAtFqtVqvt+PXq1atbt24FMGPGDE+78kIEd9lRDxQt\nLR6eop2dN2123rS05x9UezquTLmlM/WLPFU6BxzrdtoeX3YpIbtqzop1m3b9/npTo8sYrSY+PXVm\nQO8UKlpa7HDgRn1Tc+PN2MSYGR8+ofaMOvGeEYDrzfbfrG1IW1T91NK8jza+Yi7b5TIgCDJysWLF\nirq6Oq1Wu2nTpgEDQromh/R/njyRZXnlypVJSUkA8qZ8qrz5ldXDL2x+JjV7fGr2+MQpCSrPsrOy\nwgqrxZYsTfI0YPWXe5d/1lh2tvVmq6O24dq+Q9v+sOFfb99udRmWIqVZLTa51OLj+faV1WIz5Zbm\nTfkUXTKKif+bqu+rwweFz1k7c0h0hMoTdXLXjNra259+o/a/dl+Xa9uaW25VnDUXFBm/PvCFy7BA\nzCgpKWnlypVdxxw6dGjdunUAcnJyoqOj/T1FwbAhkStZlpOSkoxGoyzLuPOmMhvjcyHquSzOlEVZ\nXV25WlNUUgbgsYeGbjWOfHzyWAA1decOn/hzl2eIAyCbhS52Vostb8qnplyz1WJD54wm/tO4b/94\nEMDUt/4h/u9GqjxRdzxlBGDn/vKTVbcAvPl87Odv//OEMRkAtv/5T623bnZ+hsDLSJZlo9GYlJSk\nvK0Udrt98eLFdrs9KipqyZIlqs1VGGxI5GrBggVuHxf8/V/2ZUV66gxPW89Zytrt9gFh+FPOsMce\nGrrsuSfDwwcCqL58xmWkVhOfLKXJpdW+nW7fKN+EupLNlq8/3O+wOwaEhzXVN+/9YP+B9UeUTRfL\navd+sP/S8To/TtOV94wAHD9fAyA5YdCrc2KS4oc9N2M5gFu3W8pP73MeFrgZybKclZXV8euePXt+\n+OEHAM8++2xsbKz/JicqNiTqpKCgwGQyud0k+E4S2WxJkSZ72trWdvu+4Zq0lIiRseEA7HaHw2EH\nEHXPvV0HZ6TOtFZf891U+6issMLThwOrxdZibQFgb3d8+8cD33504MfNx5RNtRV13350oO6kmuvu\nvGcEoNHWBCA26k5d6kjnYt05l5GBm5Esyx1vMWUtA4Ann3zSPxMTHFfZUQ+Ycs3SX95ml0/cKW2W\nwzWqH2FW9op4OTjx95N+/tLPWh+IKwRws9WxKv9/7XZ7WNiA8fc/0nVwsjTJWmQrWlqs0Ym4T1/5\nz3rS1Hhz0JC/vq8dDkdba7vy86AhAyv3nr/6k9W38/PgrhkBSBihOXQShytvnai6hcHYf3i78rit\nyXWhYEBnVFJSotfr29raioqKAAwYMGDq1Kn+mprQ2JCok6qqKi9b9SmPS1GS8vPRIUdP4gyAhEGJ\n+qin/DA3L+RIuQwVXg5OADhRm32iNrv68pn121bU1J0HMCvr5YRRyV1HajXxWk2cpm64/gG9jybc\nF6Y6k5etc/9xniRJHb/W1tauXbsWwKRJk2bNmuXjqXnTnYwmPvD6jn3PtbXbH/vVldjoLfVX/7JT\nzuF6pkFAZ5SYmAigsrKysbERwOjRo51XgYcyNiTqJDPT24Xp5s+fr9frlZ83bty4fft2AOnp6e+8\n844f5uaFyWQqKCiolA97P7vlux93bPm/D2+33YoYHDn3Z7/JSJ3paWSD9dL8+e8bDIb+n2ufZWZm\netqtis4ZASgvL1caUlpamroxdSejhJFjFmV/uHHnezearXWNlgljMqovn7Xa6mKj3azOCNyMlE8M\np06dUn5NSBBrwaqKeAyJOtHr9c7lzMsmoU6YUN7hp+XDXsbs/X7z5//z7u22W7q4lDdf3uSlGynn\nvohZ6dCTjJyFh4f7bkrd0Z2Mmpqv6UYlv/ny5ncWF/7Hr/cYZhuvXa8HEDditMvIIMjozJk7C2rY\nkDqEObp8F6YQpywEcl6cCkCSpPz8fE/vMRFkZWXVyNeXGda53Xr5SpXx42y7vT3qntglv8yLGByp\nPD404h5N9AiXwQVFxpTUYfn5+b6dcR8EZUYA9n6/pbA4Nzx84O+WbBsWe9+Xu//z6wNfRA6N/vdl\nxYMGDnYeyYyCkkAfckkQkiR98803RqNR+UgrSZLRaLxw4YLg7yK9Xt9grfG0df+RHXZ7O4DrTVd/\n/+l845o5yr+v9uR1HVwpH/a+61J1QZkRgEnjHhsSEdne3va7T+auWP1z5ZTYGZkvuXQjMKMgxW9I\nFCSU83mXGda5PUTx8ZZfHe18Lovi4YnT/uWZVc6PmMt2FRQZL1y44Lw0gPqF94wUp84f/O+d7125\nehFAbPSIJ9KffyL9ly5jmFGw4qIGChIdhyjcFrvF8/7Qi2ej/uU9I8XY0Q+v+rftjdcuOxz2YbH3\n3fXZKJhwlx0Fj/z8fHPZTu87hbxrsNYUFBlFPjIR6LqTUVjYAK0mzlM3YkZBjA2JgofBYHhl6aLc\ngkW960kN1prlebMMBoOwa7eCADMiL8KNRqPacyDqN3q9/mbrjby1xtSx+sghUd3/w45Kx4/evsaM\nyBM2JAo2vah3rHR+xozILTYkCkId9a655UaY1/sdNFhr9n6/5ZMvfs1K52fMiLrism8KWitXrlQ+\nb2k1cclSWoo0WfmhwVpTWrarwVqjnO0vSZLBYFD96kehiRmRMzYkCnLK1f43bNjgcm0xpcZlZmby\nREXVMSNSsCFRqJBlueM6LixwYmJGIY4NiYiIhMDzkIiISAhsSEREJAQ2JCIiEgIbEhERCYENiYiI\nhMCGREREQmBDIiIiIbAhERGRENiQiIhICGxIREQkBDYkIiISAhsSEREJgQ2JiIiEwIZERERCYEMi\nIiIhsCEREZEQ2JCIiEgIbEhERCQENiQiIhICGxIREQmBDYmIiJwgHckAAAA8SURBVITAhkREREJg\nQyIiIiGwIRERkRDYkIiISAhsSEREJAQ2JCIiEgIbEhERCYENiYiIhMCGREREQvh/Vo24EZYveukA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "[node,elem] = squaremesh([0 1 0 1], 0.5);\n",
    "bdFlag = setboundary(node,elem,'Dirichlet');\n",
    "[elem,bdFlag] = sortelem(elem,bdFlag);\n",
    "showmesh(node,elem);\n",
    "findnode(node);\n",
    "findelem(node,elem);\n",
    "[elem2edge,edge] = dofedge(elem);\n",
    "findedge(node,edge,'all','rotvec');\n",
    "display(elem2edge);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assembling the matrix equation\n",
    "\n",
    "We discuss several issues in the assembling.\n",
    "\n",
    "### Mass matrix\n",
    "\n",
    "The mass matrix can be computed by\n",
    "\n",
    "    M = getmassmatvec(elem2edge,area,Dlambda,'RT1');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Matlab",
   "language": "matlab",
   "name": "matlab"
  },
  "language_info": {
   "codemirror_mode": "octave",
   "file_extension": ".m",
   "help_links": [
    {
     "text": "MetaKernel Magics",
     "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
    }
   ],
   "mimetype": "text/x-matlab",
   "name": "matlab",
   "version": "0.14.3"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "138px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
